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.  TESTS  OF  INITIALIZATION  PROCEDURES  WITH  THE  NRL 
^LIMITED  AREA  NUMERICAL  WEATHER  PREDICTION  MODEL 

1.  INTRODUCTION 

Various  initialization  procedures  have  been  tested  for  use  with  the  Naval 
Research  Laboratory’s  (NRL)  Limited  Area  Weather  Prediction  system.  The  system 
has  been  developed  to  study  the  development  of  eztratropical  cyclones,  which 
occurred  along  or  off  the  East  coast  of  the  U.S.  during  the  Genesis  of 
Atlantic  Lows  Experiment  (GALE)  and  the  Experiment  on  Rapidly  Intensifying 
Cyclones  over  the  Atlantic  (ERICA) .  Errors  in  the  analysis  (which  can  be  due 
to  observational  errors  and  unresolvable  scales  of  motion)  and  inaccuracies  in 
the  model  physics  give  rise  to  inertia-  gravity  wave  oscillations  in  numerical 
integrations  of  the  model.  Such  errors  are  reflected  as  unbalanced  deviations 
in  the  wind  and  mass  fields,  which  generate  freely  propagating  inertia- 
gravity  waves.  For  the  external  and  first  few  internal  vertical  modes  of  the 
numerical  model,  the  phase  speeds  of  these  free  inertia-  gravity  waves  are 
much  larger  than  the  speeds  of  meteorological  systems.  The  resulting  high 
frequency  oscillations  can  be  seen  in  the  surface  pressure  for  example,  with 
amplitudes  as  large  as  5  to  10  mb. 

Over  the  years,  various  methods  have  been  used  to  reduce  these  high 
frequency  oscillations  in  integrations  of  numerical  weather  prediction  models. 
Among  these,  initial  conditions  used  for  integrations  of  the  model  are 
sometimes  modified  or  initialized  by  application  of  various  filtering 
equations.  In  the  conventional  static  initialization  performed  on  pressure 
surfaces,  horizontal  scaling  arguments  are  used  to  derive  the  non-linear  mass 
balance  equation,  which  relates  the  geopotential  and  the  stream  function  of 
the  non-divergent  wind  for  the  larger  scale  atmospheric  motions  (see  Haltiner 
and  Williams,  1980,  for  example).  In  midlatitudes  it  has  been  customary  to  use 
the  observed  geopotential  heights  to  compute  a  statically  balanced  wind  field, 
assuming  the  wind  is  geostrophic  at  the  lateral  boundaries  of  the  model  (see 
Bengtsson,  1975,  for  example).  On  the  other  hand  in  the  tropics,  the  non- 
divergent  wind  is  used  to  derive  a  statically  balanced  temperature  from  the 
geopotential  (see  Krishnamurti,  1979,  for  example).  Since  most  numerical 
weather  prediction  models  use  vertical  coordinates  different  from  the 
pressure,  these  statically  balanced  mass  and  wind  fields  must  be  then 
Manuscript  approved  May  2.  1990 
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Interpolated  to  the  model  vertical  coordinates.  Introducing  some  noise  in  the 
initial  conditions.  Sundqvist  (1975)  used  horizontal  scaling  arguments  to 
derive  the  non-linear  mass  balance  equation  directly  on  the  vertical  sigma 
levels  of  his  numerical  model.  In  the  normal  mode  initialization,  used  in 
global  numerical  weather  prediction  models  (Andersen,  1977,  Daley,  1979, 
Temperton  and  Williamson,  1981,  Williamson  and  Temperton,  1981,  for  example), 
the  normal  modes  of  the  numerical  model  are  computed  and  by  projecting  the 
initial  wind  and  mass  fields  (which  have  been  interpolated  to  the  model 
coordinates)  onto  these  normal  modes,  the  high  frequency  inertia-  gravity 
waves  can  then  be  removed.  However  in  limited  area  models,  it  is  not  possible 
to  define  the  horizontal  structure  of  the  normal  modes.  In  the  vertical  mode 
initialization  scheme  of  Bourke  and  McGregor  (1983),  filtering  conditions  for 
the  inertia-  gravity  waves  are  applied  to  the  model  dynamical  equations  to 
derive  linear  diagnostic  equations  for  the  mass  divergence  and  generalized 
geopotential  for  the  first  few  vertical  modes  of  the  model.  With  the  further 
condition  that  the  linearized  potential  vorticity  is  unchanged  by  the 
procedure,  the  equations  can  be  solved  iteratively  for  the  amplitude  of  the 
high  frequency  gravity  waves  in  the  initial  conditions.  The  scheme  has  been 
shown  to  be  an  application  of  the  normal  mode  initialization  scheme  used  in 
global  models,  without  the  horizontal  structure  of  the  normal  modes  having  to 
be  computed  (Juvanon  du  Vachat,  1986;  Temperton,  1988). 

To  reduce  the  amplitude  of  these  high  frequency  gravity  wave  oscillations 
in  the  NRL  model,  several  basic  types  of  initialization  procedures  have  been 
used.  A  static  initialization  procedure,  which  had  been  developed  for  the  NRL 
model,  is  tested.  With  the  future  implementation  of  a  wind  profiling  network 
in  the  U.S.,  there  is  renewed  interest  in  deriving  geopotential  heights  from 
the  wind  field  in  the  midlatitudes.  In  the  NRL  scheme  then,  the  p  ,n-divergent 
wind  is  first  computed  from  the  analyzed  winds  on  the  pressure  surfaces.  The 
computed  non-divergent  wind  and  the  analyzed  temperatures  are  then 
interpolated  to  the  model  vertical  coordinates.  A  diagnostic  relation  for  the 
geopotential  on  the  sigma  surfaces  of  the  numerical  model  is  derived,  by 
setting  the  tendency  of  mass  divergence  to  zero  and  ignoring  vertical 
advection  and  friction  in  the  divergence  equation.  For  the  boundary  condition 
on  the  geopotential,  we  generalized  the  conventional  geostrophic  relationship. 
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deriving  diagnostic  conditions  for  the  normal  derivatives  of  the  geopotential 
by  ignoring  the  tendencies  of  momentum  in  the  momentum  equations. 

A  second  type  of  initialization  procedure,  based  on  the  vertical  mode 
initialization  scheme  of  Bourke  and  McGregor  (1983),  has  also  been  developed 
for  the  NHL  model.  As  is  customary,  we  keep  the  geopotential,  temperature  and 
pressure  fixed  at  the  lateral  boundaries  of  the  model.  To  provide  a  boundary 
condition  for  the  mass  divergence  however,  an  approximate  mass  divergence  is 
computed  along  the  lateral  boundary  using  the  thermodynamic  equation.  In  our 
scheme,  changes  in  the  tangential  wind  along  the  lateral  boundaries  are 
consistent  with  the  changes  in  the  vorticity  and  divergence  computed  by  the 
scheme. 

The  effect  of  the  split-explicit  scheme,  which  is  used  for  integration  in 
time  in  the  NRL  model,  and  the  non-linear  initialization  procedures  in 
reducing  gravity  wave  oscillations  in  integrations  of  the  NRL  model  are 
investigated.  The  influence  of  two  different  lateral  boundary  treatments  and 
two  different  grids  of  differing  resolution  and  domain  size  are  investigated. 
In  these  integrations,  idealized  boundary  conditions  are  obtained  by 
interpolation  from  operational  analyses.  To  minimize  the  impact  of  noise  from 
the  boundaries  influencing  the  interior,  the  boundary  conditions  are  derived 
from  initialized  fields  for  the  cases  of  integrations  starting  from 
initialized  fields. 

The  numerical  model,  the  vertical  modes,  the  split-explicit  integration 
scheme  and  the  two  different  lateral  boundary  treatments  are  described  in 
section  2.  In  section  3,  the  static  initialization  procedure  is  described  and 
illustrated.  The  vertical  mode  scheme  is  then  described  in  section  3.  The 
convergence  of  the  scheme  is  also  shown  for  the  two  different  grids.  A  low 
resolution  grid  with  a  resolution  of  2°  in  longitude  by  1.5°  in  latitude, 
covers  the  continental  U.S.,  and  a  high  resolution  grid  of  0.5°  resolution  in 
latitude  and  longitude  covers  the  eastern  U.S.  In  section  5,  integrations  with 
the  split-explicit  scheme  for  time  integration  are  compared  with  a  centered 
difference  scheme  on  the  high  resolution  grid,  for  uninitialized  initial 
conditions.  Integrations  with  the  split-explicit  scheme  are  then  compared  for 
varying  degrees  of  static  initialization.  In  section  6,  integrations  with 
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initial  conditions,  initialized  with  the  vertical  mode  scheme,  are  compared  to 
integrations  with  uninitialized  fields  for  both  the  low  and  high  resolution 
grids.  Differences  due  to  the  different  lateral  boundary  treatments  are 
compared.  The  results  are  summarized  and  discussed  in  section  7. 
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2. 


MODEL  DESCRIPTION 


The  Naval  Research  Laboratory’s  primitive  equations  model  (Madala  et  al., 
1987)  uses  sigma  coordinates  in  the  vertical  and  incorporates  topography  and 
physical  parameterizations  of  the  boundary  layer  and  precipitation  processes. 
The  model  is  integrated  in  time  using  the  efficient  split-explicit  method 
(Madala,  1981).  In  the  horizontal,  an  Arakawa  C  grid  (Arakawa  and  Lamb,  1977) 
is  used  with  a  latitude  and  longitude  grid.  The  finite  difference  scheme  is  a 
second  order  quadratic  conserving  scheme.  The  model  has  been  used  for  example 
to  study  the  east  coast  snow  storm  of  10-12  February,  1983  (Chang  et  al., 
1989).  For  experiments  in  this  paper,  ten  layers  of  equal  thickness  (with 
A a  ■  0.1)  are  used  in  the  vertical  sigma  ( a )  coordinate  from  the  surface 
(<r  «=  1)  to  the  model  top  (a  *  0).  The  model  includes  large  scale  precipitation 
and  a  cumulus  parameterization  using  a  modified  Kuo  scheme.  Unstable  lapse 
rates  are  removed  by  a  dry  convective  adjustment  scheme  following  Manabe  et 
al.  (1965).  The  boundary  layer  is  parameterized  using  a  drag  coefficient 
formulation,  and  second  order  horizontal  diffusion  is  included. 

Two  different  model  grids  are  used  for  the  model.  A  low  resolution  grid 
(called  the  US  grid),  with  a  resolution  of  2°  longitude  by  1.5°  latitude, 
covers  the  continental  U.S.  in  a  domain  from  140. 0°W  to  40.0®W  and  10.0°N  to 
70 . 0°N .  The  other  grid  (called  the  GALE  grid)  is  a  high  resolution  grid  of 
0.5°  resolution  in  latitude  and  longitude,  covering  the  eastern  U.S.  from 
102. 5®W  to  57.5®W  and  22.5°N  to  47.5®N.  Analyses  at  14  standard  pressure 
levels  (from  50  to  1000  mb)  on  a  2.5°  hemispheric  grid  from  the  National 
Meteorological  Center  (NMC)  provide  the  initial  conditions  and  idealized 
boundary  conditions  for  model  integrations  in  this  paper.  The  NMC  2.5® 
resolution  hemispheric  analysis  is  interpolated  to  the  horizontal  model  grid 
using  Lagrangian  cubic  polynomial  interpolation.  The  thermodynamic  variables 
are  interpolated  in  the  vertical  to  the  model  sigma  levels  assuming  they  are 
linear  in  the  log  of  the  pressure,  while  the  wind  components  are  interpolated 
assuming  they  are  linear  in  the  pressure.  An  enveloped  topography  is  derived 
for  each  model  grid  from  the  U.S.  Navy’s  global  10  minute  elevation  data,  by 
computing  the  average  height  for  each  model  grid  square  and  adding  one 
standard  deviation.  On  each  model  grid,  the  enveloped  topography  is  then 
filtered  by  using  the  two-dimensional  triangular  smoother-desmoother  of 
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Shapiro  (1970),  to  filter  out  any  wavelengths  in  the  topography  of  twice  the 
grid  distance.  For  use  with  the  vertical  mode  initialization  scheme,  the  model 
topography  is  further  smoothed  by  using  one  pass  of  the  two-dimensional  nine 
point  triangular  smoother  (Shapiro,  1970).  Climatalogical  mean  sea  surface 
temperatures  for  the  month  of  January  of  one  degree  resolution,  taken  from 
Reynolds  (1982),  are  interpolated  to  the  model  grids. 


2.1  Vertical  Modes  of  the  Model 


In  order  to  solve  for  the  vertical  modes  of  our  numerical  model,  we 
linearize  about  a  basic  state  at  rest  with  a  mean  temperature  T*.  which  is  a 
function  of  sigma  only,  and  separate  the  linear  and  non-linear  terms  in  the 
model  dynamical  equations.  The  model  dynamical  equations  in  flux  form  can  then 
be  written  in  matrix  notation  as 


0P3U 

dt 

8t 

8p£ 

8t 

8p 

*s 

8t 


V 


n2td 


0 


(2.1) 

(2.2) 

(2.3) 

(2  .U) 


and  the  hydrostatic  relation  and  continuity  equation  are  written 

*  -  *e  -  ‘2-5> 

P8i  -  D  (2.6) 

where  column  vectors  are  in  bold  type,  Mj_,  M2,  are  matrices,  and  N2T 

represents  the  transpose  of  vector  N2.  The  vertical  sigma  coordinate  (a)  is 
defined  by  the  ratio  of  the  pressure  p  to  the  surface  pressure  ps  (Phillips, 
1957).  The  dynamical  variables  are  in  vector  form,  where  the  elements  of  the 
vectors  represent  the  values  of  the  variables  at  each  of  the  ten  model  sigma 
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levels  for  a  single  horizontal  grid  point.  The  vectors  u,  ▼  represent  the 
horizontal  wind,  T  the  temperature  and  ^  the  geopotential,  at  each  of  the 
model  sigma  levels,  which  are  defined  at  the  middle  of  each  of  the  sigma 
layers  in  the  vertical.  is  the  surface  geopotential  (at  a  -  1).  The 

generalized  geopotential  •  is  defined  as 

♦  *  Pst^s  +  RT*-**]  (2.7) 

where  the  average  geopotential  f*  on  the  sigma  surface  is  related  to  the  mean 
temperature  T*  through  the  hydrostatic  relation  T*.  The  vertical 

motion  a  in  the  sigma  coordinate  (a)  is  staggered  in  the  vertical,  being 
defined  at  sigma  levels  at  the  boundaries  between  the  vertical  layers. 
Subscripts  representing  the  horizontal  grid  points  on  the  C  grid  (see  also 
grid  mesh  in  Fig.  2  in  section  2.3)  have  been  omitted  for  clarity.  The  x  and  y 
coordinates  are  defined  by  multiplying  the  longitude  and  latitude  in  radians 
by  the  average  radius  of  the  earth.  The  mass  divergence  D  on  the  sigma 
surfaces  is  defined  on  the  C  grid  in  our  spherical  coordinates  by 


i  V 

y 


h  p' 
y 


u  ) 


i  6  (  V 


(2.8) 


Here  the  difference  operator  6  is  defined  in  the  x  direction,  using  the 
generalized  geopotential  •  as  an  example,  by 


*(x+Ax/2)  -  *(x-Ax/2) 

h  Ax 
x 


(2.9) 


where  Ax  is  the  grid  spacing  for  the  x  coordinate  and  hx  (equal  to  the  cosine 
of  the  latitude  for  our  coordinate)  is  the  map  factor  for  the  x  coordinate.  A 
similar  difference  operator  is  defined  for  the  y  coordinate,  where  the  map 
factor  is  hy  ■  1  in  our  case.  An  averaging  operator  is  also  defined  in  the  x- 
direction,  using  the  surface  pressure  as  an  example,  by 

p  (x+Ax/2)  +  p  (x-Ax/2) 

p6  -  - $ -  (2-10) 


A  similar  averaging  operator  is  defined  for  the  y-coordinate  and  a  two 
dimensional  averaging  can  be  defined  as 
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(2.11) 


Elements  of  the  matrices  M^,  M2,  and  the  vector  N2  are  functions  of  sigma 
only.  The  vectors  on  the  right  hand  side  (RHS)  of  Eqs.  (2.1),  (2.2)  and  (2.3) 
include  Coriolis,  friction,  non-linear  advection  and  diabatic  terms.  Details 
of  the  vector  and  matrix  elements  can  be  found  in  the  report  by  Madala  et  al. 
(1987) . 


Solutions  to  the  homogeneous  equations  (in  which  terms  on  the  RHS  of  Eqs. 
(2.1),  (2.2),  (2.3)  and  (2.4)  are  zero)  are  freely  propagating  gravity  waves . 
Now  eliminating  all  variables  in  the  homogeneous  equations  except  ♦,  we  obtain 


Mj  V2*  -  0 


(2.12) 


where  the  matrix  M3  is  defined  as  M3  =  M^  M2  +  (RT*-^*)N2T  and  whose  elements 
are  only  functions  of  the  vertical  sigma  coordinate  and  do  not  depend  on  the  x 
and  y  coordinates.  In  our  x  and  y  spherical  coordinates,  the  two  dimensional 
Laplacian  V2  is  defined  by 


V2  * 


i  v 

y 


h  6  t  ) 
y  x 


r  Si  h l  6_ 


*  ) 


(2.13) 


Similarly,  we  can  show  that  the  mass  divergence  D  satisfies  the  same  equation 
(2.12),  while  the  vorticity  of  the  background  flow  is  not  affected  by  gravity 
wave  motions.  The  equation  (2.12)  is  separable  and  by  separating  the  vertical 
structure,  a  set  of  vertical  eigenmodes  and  corresponding  eigenvalues  can  be 
obtained  (see  also  Kasahara  and  Puri,  1981,  for  example).  In  terms  of  our 
matrix  notation,  the  eigenvectors  and  corresponding  eigenvalues  are 

found  by  solving  the  matrix  equation 

M3  «k  “  *k  *k-  (2.14) 

If  E  represents  the  eigenvector  matrix  (with  each  column  representing  an 
eigenvector  e^)  of  matrix  M3,  and  A  is  the  diagonal  matrix  with  the  diagonal 
elements  given  by  the  eigenvalues  ,  then  we  can  write 

M3  E  -  A  E.  (2.15) 
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and  by  multiplying  Eq.  (2.12)  by  the  inverse  of  E  we  have 


a2 

0  e 


8t' 


a  r 


(2.16) 


where  e  =  E"1  •  ,  and  we  have  used  the  property  of  the  eigenvector  matrix  that 
E-1  M3  E  =  A  ,  which  can  be  easily  verified  by  comparing  the  elements  of 
matrices  M3  E  and  E  A  (see  Strang,  1988,  p254,  for  example).  For  each 

vertical  mode  we  have 


a2 

8  e. 


8t' 


X.  V2  e. 

k  k 


(2.17) 


where  e^  is  the  amplitude  of  the  generalized  geopotential  and  X^  is  the 


TABLE  1:  Mean  temperatures  for  12Z  January  23,  1986,  on 
sigma  surfaces  of  10  layer  model  for  GALE  grid. 


Model  Level 

Sigma  Level 

Temperature 

(°K) 

1 

0.05 

205.9 

2 

0.15 

214.0 

3 

0.25 

222.1 

4 

0.35 

236.1 

5 

0.45 

249.4 

6 

0.55 

259.4 

7 

0.65 

266.6 

8 

0.75 

271.9 

9 

0.85 

275.9 

10 

0.95 

278.4 

9 


TABLE  2:  The  equivalent  depths  and  the  phase  speeds  for 
the  vertical  modes  of  a  ten  layer  model  for  the 
mean  temperature  profile  defined  in  Table  1. 


Mode  No.  Equiv.  Depth  Phase  Speed 

(meters)  (meters  s-*) 


1 

9,399.0 

303.5 

2 

1,508.0 

121.6 

3 

226.6 

47.1 

4 

85.7 

29.0 

5 

30.3 

17.2 

6 

13.7 

11.6 

7 

6.4 

7.9 

8 

2.7 

5.1 

9 

0.9 

3.0 

10 

0.2 

1.3 

eigenvalue  for  the  kth  mode.  Eq.  (2.17)  is  a  wave  equation  describing  the 
propagation  of  the  free  gravity  mode  whose  phase  speed  c^  is  given  by  fX^. 
In  the  linearization  used  here  to  obtain  the  gravity  modes,  the  Coriolis  term 
is  combined  with  the  non-linear  terms,  and  as  such  the  gravity  modes  are 
applicable  only  for  high  frequencies  or  small  time  periods  compared  to  the 
period  of  inertial  oscillations.  However  it  can  be  noted  that  the  vertical 
structure  of  the  modes  derived  in  this  case  is  unchanged  even  if  the  Coriolis 
terms  are  included  with  the  linear  terms  on  the  left  hand  side  of  Eqs.  (2.1) 
and  (2.2).  In  this  case  the  linearization  defines  the  inertia-  gravity  modes 
(see  section  4.1) . 

As  an  example  the  eigenmodes  are  computed  for  the  ten  layer  model  using  a 
emperature  profile  taken  from  the  NMC  2.5°  hemispheric  analysis  for  12Z 
January  23,  1986  during  the  second  Intensive  Observing  Period  (IOP)  of  GALE. 
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The  temperatures  are  Interpolated  to  the  model  grid  and  averaged  on  the  sigma 
surfaces  over  the  GALE  grid  covering  the  eastern  U.S.  The  average  temperatures 
on  the  model  sigma  levels  at  the  middle  of  each  layer  are  shown  in  Table  1. 
The  eigenvalues  and  the  phase  speeds  for  each  of  the  ten  vertical  modes 
computed  are  shown  in  Table  2.  The  eigenvalues  are  given  in  terms  of 
equivalent  depths  defined  by  X^/g  .  where  g  is  the  acceleration  due  to 
gravity.  The  corresponding  vertical  structure  for  each  of  the  modes  is  shown 
in  Fig.  1.  The  number  of  zero  crossings  in  the  vertical  structure  of  each  mode 
is  given  by  the  mode  number  minus  one.  The  phase  speeds  of  the  modes  vary 
from  303.5  m  s"1  for  the  external  (first)  mode  to  1.3  m  s"1  for  the  tenth 
mode.  For  the  first  three  modes,  the  external  mode  and  first  two  internal 
modes,  the  phase  speeds  are  much  faster  than  typical  meteorological  systems, 
which  are  typically  less  than  20  m  s”1.  Then  to  integrate  a  model  with  a 
centered  difference  scheme  in  time  and  a  grid  spacing  Ax,  the  time  step  2  At 
must  be  small  enough  to  satisfy  the  CFL  condition  2  At  (U+c)  /Ax  £  1  ,  for 
the  external  gravity  mode  which  has  the  fastest  phase  speed  c,  and  where  U  is 
a  maximum  advecting  wind  speed  (see  for  example  Mesinger  and  Arakawa,  1976, 
p51) .  To  allow  the  NHL  model  to  be  integrated  more  efficiently  at  a  larger 
time  step,  more  appropriate  for  meteorological  systems,  the  split-explicit 
scheme  was  developed. 

2.2  Split-Explicit  Time  Integration 

In  the  split-explicit  method  of  Madala  (1981),  a  centered  difference 
scheme  is  used  with  a  large  time  step  to  compute  initial  estimates  of  the 
tendencies  of  the  model  variables  for  all  the  terms,  except  for  diffusion 
which  uses  forward  differencing.  The  time  step  satisfies  the  CFL  condition  for 
the  phase  speed  of  the  4^  gravity  mode.  To  step  the  model  variables  at  time 
t  -  At  forward  in  time  by  2  At,  using  the  centered  difference  scheme,  the  non¬ 
linear  terms  and  forces  are  computed  at  time  t.  These  first  estimates  of  the 
tendencies  of  the  mass  (surface  pressure  ps,  temperature  T)  and  wind  fields 
(u,  v)  are  then  corrected  for  the  motion  of  the  higher  frequency  gravity 
modes,  assuming  that  the  computed  non-linear,  Coriolis,  diabatic  and  friction 
forcing  terms  are  constant  during  the  time  step  of  2  At.  The  specific  humidity 
q  is  not  corrected.  To  obtain  the  corrections,  the  amplitudes  of  the 
deviations  of  the  divergence  D  -  D(t)  and  the  generalized  geopotential 
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♦  -  ♦ ( t )  are  Integrated  over  the  interval  of  2  At  at  smaller  time  steps,  for 
each  of  the  first  three  vertical  modes.  For  the  external  mode  a  time  step  of 
At/8  is  used,  while  for  the  first  two  internal  vertical  modes  time  steps  of 
At/4  and  At/2,  respectively,  are  used.  The  average  of  these  deviations,  over 
the  interval  of  twice  the  large  time  step  2  At,  is  then  used  to  correct  the 
initial  explicit  estimate  of  the  variables.  Further  details  can  be  found  in 
the  Appendix.  For  the  integration  of  the  deviations  of  the  divergence  a 
lateral  boundary  condition  is  also  required  for  the  generalized  geopotential. 
A  boundary  value  for  •  -  • ( t )  is  computed  by  a  linear  interpolation  from  the 
boundary  values  at  t-At  and  t.  Further  pragmatic  boundary  conditions  are 
provided  by  reducing  the  amplitude  and  phase  of  the  deviations  of  the 
divergence  and  generalized  geopotential  in  a  boundary  zone  (see  the  Appendix). 
Besides  providing  a  4  to  5  times  saving  in  computer  time  over  explicit 
methods,  the  averaging  of  the  lower  gravity  wave  eigenmodes  can  be  expected  to 
reduce  the  amplitude  of  the  freely  propagating  higher  frequency  gravity  waves. 

2.3  Lateral  Boundary  Conditions 

To  update  the  large  scale  flow  at  the  horizontal  boundaries  during  the 
integration  of  the  NRL  model,  externally  prescribed  boundary  conditions  are 
required  for  both  the  u  and  v  components  of  the  wind  field  and  the  mass  and 
humidity  fields.  In  the  model,  the  mass  and  humidity  variables  of  surface 
pre' sure  ps,  temperature  T,  geopotential  f  (or  generalized  geopotential  •), 
and  specific  humidity  q  are  defined  at  the  lateral  boundary.  The  wind 
components  are  staggered  in  the  C  grid  and  as  applied  to  the  NRL  model,  the 
tangential  wind  is  defined  at  the  boundary,  while  the  normal  wind  is  staggered 
half  a  grid  point  in  from  the  boundary  (see  Fig.  2).  In  this  paper  two 
different  lateral  boundary  treatments,  the  tendency  relaxation  scheme  of 
Perkey  and  Kreitzberg  (1976)  and  the  Davies  (1976,  1983)  relaxation  scheme, 
are  used  and  compared  for  use  with  various  initialization  procedures.  To 
provide  the  model  boundary  conditions  in  our  experiments,  idealized  boundary 
values  and  tendencies  are  derived  from  the  NMC  2.5  degree  hemispheric 
analyses,  interpolated  to  the  model  grid. 

(a)  Perkey  Kreitzberg  scheme. 

In  the  Perkey  Kreitzberg  scheme,  model  computed  tendencies  for  each  of  the 
dependent  variables  are  blended  with  specified  boundary  tendencies  in  a 
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boundary  zone  of  5  points.  After  each  time  step,  the  model  tendencies  for  each 
of  the  independent  model  variables  are  adjusted  according  to 


8a 

8a  ‘ 

8a 

—  -=  (l-a) 

— 

+  a 

— 

8t 

.  8t  . 

m 

.  8t  . 

where  a  is  the  independent  model  variable  psu,  Psv»  Ps*  Ps*»  or  Ps<l» 
subscript  m  represents  the  model  computed  value,  and  subscript  b  the 
prescribed  boundary  value,  a  is  a  linear  function  of  the  minimum  distance  (n) 
from  the  lateral  boundary,  in  units  of  the  grid  spacing.  As  in  Chang  et  al. 
(1989),  we  use  for  the  mass  and  humidity  variables  (that  is,  for  p8,  psT,  psq) 


1  -  —  for  n  £  5 
0  for  n  >  5 


(2.19) 


■■  y  y 

while  for  the  wind  components  psu,  p8v,  which  are  staggered  on  the  grid  mesh, 
we  use 


a 


u 


a 


V 


I 


1 


l 

5y 


for  n  £  0.5 
for  n  l  1 
for  n  £  0.5 
for  n  1  1 


(2.20) 


(2.21) 


The  boundary  tendencies,  derived  from  the  twelve  hourly  NMC  analyses,  are 
updated  every  twelve  hours  of  model  integration. 

(b)  Davies  scheme. 

In  the  Davies  scheme,  computed  model  variables  are  relaxed  to  the  boundary 
values  themselves  at  each  time  step  in  a  boundary  zone  of  6  points.  In  this 
case,  the  model  variables  are  adjusted  at  each  time  step  according  to 

a  -  (l-a)  +  a  ab  (2.22) 

where,  following  Gr^nas  et  al.  (1987),  we  define  a  as  a  quadratic  function  of 
the  minimum  distance  (n)  from  the  lateral  boundary,  in  units  of  the  grid 
spacing.  For  the  mass  variables  we  use 
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1 

for 

n  ■ 

0 

a 

for 

1  i 

n  £  5 

(2.23) 

0 

for 

n  2 

6 

The  functions  au  and  av  for 

the  u 

and 

v  components  of 

the  wind  field  are 

defined  similarly  as  in  Eqs. 

(2.20) 

and 

(2.21),  using  the 

averaging  operators 

in  the  x  and  y  coordinate  directions,  respectively.  At  each  time  step,  the 
boundary  values  are  computed  by  a  linear  interpolation  in  time  from  the  12 
hourly  NMC  analyses. 
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3. 


STATIC  INITIALIZATION 


In  the  static  initialization  procedure  used  at  NRL,  a  diagnostic 
relationship  is  derived  for  the  geopotential  on  the  sigma  surfaces  of  the 
model.  The  non-divergent  wind  is  first  computed  for  the  analyzed  winds  on  the 
pressure  surfaces.  The  NRL  model  differs  from  many  other  limited  area  models 
in  that  the  tangential  wind  is  specified  at  the  lateral  boundary  of  the  model, 
instead  of  the  normal  wind.  To  compute  the  non-divergent  wind  then  we  solve 
for  the  stream  function,  assuming  that  the  tangential  wind  along  the  lateral 
boundary  is  purely  non-divergent.  The  non-divergent  wind  and  the  analyzed 
temperatures  are  then  interpolated  to  the  sigma  coordinates  of  the  model.  A 
diagnostic  relation  is  then  derived  for  the  geopotential  on  the  sigma  surfaces 
of  the  numerical  model,  by  ignoring  the  tendency  of  divergence,  vertical 
advection  and  friction.  The  initial  non-divergent  wind  and  analyzed 
temperatures,  interpolated  to  the  sigma  surfaces,  are  used  to  compute  the  non¬ 
linear  forcing  terms.  Boundary  conditions  for  the  normal  derivatives  of  the 
geopotential  are  obtained  by  ignoring  the  tendencies  in  the  momentum 
equations . 

3.1  Non-Divereent  Wind  on  Pressure  Surfaces 


For  large  scale  atmospheric  motions,  the  divergence  of  the  velocity  field 
is  an  order  of  magnitude  smaller  than  the  vorticity.  To  a  first  approximation 
then,  the  flow  can  be  considered  non-divergent  on  surfaces  of  constant 
pressure.  On  pressure  surfaces,  the  vorticity  Jp  and  divergence  Dp  are  defined 
on  the  model  horizontal  grid  as 


h  *„< 


h  v 
7 


1 


D 

P 


h  *x' 
7 


h  u 
7 


r  6  1  h?v 


(3.1) 


(3.2) 


where  u  and  v  are  the  analyzed  wind  components  on  a  pressure  surface, 
interpolated  to  the  horizontal  grid  of  the  model.  The  non-divergent  flow  can 
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then  be  described  by  introducing  a  stream  function  f,  so  that  the  non- 
divergent  wind  components  u^  and  at  that  pressure  level  are  given  by 

u^  -  -  6yf  (3.3) 

vj»  - 

By  computing  the  vorticity  on  the  pressure  surface,  we  can  solve  Poisson’s 
equation, 

V2  f  -  £p  (3.4) 

for  the  stream  function  f,  and  whence  for  the  non-divergent  wind.  To  provide 
boundary  conditions  for  the  stream  function,  we  assume  that  the  analyzed 


tangential 

wind 

along  the 

boundary  of 

our  model  domain 

is  purely  non- 

divergent. 

Then 

we  obtain 

the  Neumann 

boundary  conditions 

for  the  stream 

function 

6yf  “  -u 

at 

y  -  yo-  71 

(3.5) 

1 

< 

at 

*  -  *0*  *1 

where  jq,  y^  give  the  southern  and  northern  lateral  boundaries  and  xq.  x^  give 
the  western  and  eastern  boundaries.  In  this  case  the  stream  function  is  not 
unique,  since  adding  a  constant  value  to  the  solution  of  the  stream  function 
is  also  a  solution.  As  described  in  the  report  of  Sashegyi  and  Madala  (1989), 
to  obtain  a  unique  solution  for  the  stream  function,  we  prescribe  a  value  for 
the  stream  function  f  *  0  at  a  single  arbitrary  point  on  the  boundary  and  use 
the  elliptic  equations  solver  of  Madala  (1978,  1981).  Choice  of  a  zero  value 

as  a  first  guess  of  the  solution  of  the  stream  function  then  leads  to  the 
efficient  convergence  of  Madala' s  elliptic  equations  solver. 

As  eloquently  described  by  Lynch  (1989),  the  partitioning  of  the  wind  into 
the  divergent  and  non-divergent  parts  is  not  unique  in  a  limited  area  domain. 
The  choice  of  using  the  tangential  wind  on  the  boundary  to  describe  the  normal 
gradient  of  the  stream  function,  minimizes  the  kinetic  energy  in  the  divergent 
wind  and  does  not  lead  to  a  mixed  divergent  and  non-divergent  term  in  the 
kinetic  energy  balance.  It  can  be  noted  that  choice  of  defining  the  normal 
wind  at  the  boundary  in  some  models,  instead  of  the  tangential  wind  used  in 
the  NRL  model,  leads  to  solving  Poisson’s  equation 
V2  X  m  Dp  (3.6) 

for  a  velocity  potential  f,  which  defines  the  divergent  component  of  the  wind. 
The  non-divergent  wind  is  then  computed  by  subtracting  the  divergent  component 
from  the  original  wind  field.  For  the  boundary  conditions  in  this  case,  it  is 
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assumed  that  x  m  0  at  the  lateral  boundaries,  that  is,  It  is  assumed  that  the 
divergent  component  of  the  wind  is  zero  on  the  boundary. 


3.2  Static  Non-Linear  Mass  Balance 

For  large  scale  atmospheric  motions,  the  time  tendency  of  the  divergence 
is  small  compared  to  the  other  terms  in  the  divergence  equation.  A  diagnostic 
equation  for  the  geopotential  on  the  sigma  surfaces  of  the  model  can  then  be 
derived  by  ignoring  this  term.  Ve  can  rewrite  the  model  momentum  equations 
(2.1)  and  (2.2)  in  the  form 
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Bp^v 

*s 

+ 

6  4  + 

R  Fff  pe 

N 

(3.8) 

8t 

8 

yT 

y  s 

V 

where  we 

have 

included 

the  non-linear 

advection 

of  momentum,  the  Coriolis 

force  and  friction  in  the  terms  on  the  RHS.  The  equation  for  the  mass 
divergence  on  the  sigma  surfaces  is  then  obtained  by  taking  6X  of  hy  times  Eq. 
(3.7)  and  adding  6y  of  hx  times  Eq.  (3.8)  giving 


SD 

—  +  V*[p8V(^8)]  -  nd  -  7*[RT  Vp8  +  p8V^s]  (3.9) 


where  the  horizontal  geopotential  gradient  is  given  by 
we  have  utilized  the  two  dimensional  divergence  operator  V 


V* 


R-  rS(hN)+  r^Ch7**) 
v  n  x  y  u  n  v  x  v 

7  7  * 


vmpsv#>  -  J  «,<  hy  V,  ej  >  «•  i  «yc  ^  V  ’ 


*  (^x^*  ^y^)  and 

to  define 


(3.10) 

(3.11) 


Here  we  have  also  defined  the  vector  Ny  by  Ny  *  (Nu,  Hy).  Since  the  time 
tendency  of  the  divergence  is  assumed  small  compared  to  the  other  terms  in  the 
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mass  divergence  equation  (3.9),  we  can  ignore  this  term.  That  is,  setting  the 
tendency  of  the  mass  divergence  zero  in  Eq.  (3.9),  we  have 


V*[P8V(#-*S)]  -  Nd  -  V*[RT  Vps  +  psV*s] 


(3.12) 


To  compute  the  term  Nj),  we  first  compute  the  non-divergent  wind  on  the 
pressure  surfaces,  as  described  in  section  3.1.  The  non-divergent  wind  and  the 
analyzed  temperatures  are  then  interpolated  to  the  sigma  surfaces  of  the  model 
(see  the  beginning  of  section  2  on  the  model  description).  The  non-divergent 
wind  is  then  used  to  compute  the  terms  on  the  RHS  of  Eqs.  (3.7)  and  (3.8) 
ignoring  vertical  advection  and  friction  (leaving  horizontal  advection  and 
Coriolis  forces).  The  analyzed  temperatures  and  a  surface  pressure  computed  on 
the  model  topography  by  interpolation  are  used  to  compute  the  remaining 
forcing  terms  on  the  RHS  of  Eq.  (3.12). 


To  solve  Eq.  (3.12),  boundary  conditions  are  required  for  the 
geopotential.  To  obtain  the  boundary  conditions  we  ignore  the  tendencies  of 
the  u  and  v  momentum  in  Eqs.  (3.7)  and  (3.8)  to  obtain  the  Neumann  boundary 
conditions  for  the  geopotential 


P*  6  (4-4  )  -  N  -  R  T x6  p  -  p*  6  4 
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where  y  ■  yg,  y^  give  the  southern  and  northern  lateral  boundaries  and  x  *  xq, 
xi  give  the  western  and  eastern  boundaries.  In  Fig.  2,  our  boundary  conditions 
for  the  geopotential  gradient  are  defined  along  the  lines  given  by  the  shorter 
dashes,  half  a  grid  point  in  from  the  lateral  boundary.  Essentially,  we  have 
generalized  the  conventional  geostrophic  boundary  condition  relating  the 
gradient  of  the  geopotential  to  the  geostrophic  wind,  by  using  the  non- 
divergent  wind  and  including  the  non-linear  horizontal  advection  term.  The 
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elliptic  equation  (3.12)  can  now  be  solved  for  the  geopotential  using  the 
Neumann  boundary  conditions  (3.13)  and  (3.14)  with  an  elliptic  equations 
solver,  such  as  the  Stabilized  Error  Vector  Propagation  (SEVP)  solver  (see 
Madala,  1978;  Sashegyi  and  Madala,  1989).  For  the  SEVP  solver,  the  initial 
geopotential  derived  from  the  observed  temperatures  is  used  to  provide  the 
prescribed  value  at  a  point  on  the  lateral  boundary  (to  give  a  unique 
solution)  and  to  provide  the  first  guess  for  the  solver. 


The  balanced  temperatures  on  the  sigma  surfaces  are  then  found  from  the 
balanced  geopotential,  by  inverting  the  hydrostatic  equation.  Use  of  the 
hydrostatic  equation  to  compute  the  temperature  introduces  a  2  A a  saw  tooth 
wave  structure  in  the  vertical  temperature  profile.  To  remove  this  buckling, 
deviations  of  the  temperatures  t'  from  the  mean  temperature  T*  are  computed  at 
sigma  levels  at  the  boundaries  of  the  vertical  layers  (half  way  between  the 
model  sigma  levels). 


T  - 


T*  - 


(3.15) 


The  corrected  temperatures  are  then  obtained  by  adding  the  deviations 
(interpolated  back  to  the  model  sigma  levels)  to  the  mean  temperature  and 
removing  the  mean  of  the  deviations,  so  that  the  mean  is  unchanged. 

T  -  T*  +  t’  -  <  t’  >  (3.16) 


where  <  >  is  the  horizontal  average  over  the  model  domain  on  the  sigma 

surface. 


3 .3  Application  of  Static  Initialization 

As  an  illustration  of  the  scheme,  the  analyzed  winds  and  temperatures, 
taken  from  the  NMC  analysis  for  12Z  January  23,  1986,  are  first  interpolated 
to  the  model  horizontal  GALE  grid  for  each  of  the  14  standard  pressure  levels 
from  50  to  1000  mb  in  a  domain  covering  the  eastern  U.S.  The  horizontal 
resolution  of  the  GALE  grid  is  0.5  in  latitude  and  longitude  in  a  domain 
102. 5°W  to  57.5°W  and  22.5°N  to  47.5°N.  The  non-divergent  wind  is  then 
computed  as  outlined  above  in  section  3.1.  In  Fig.  3a  we  show  the  root-mean- 
square  (rms)  changes  that  result  in  the  u-component  of  the  wind  field  after 
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computing  the  non-divergent  wind.  Changes  in  the  v-component  of  the  wind  field 
are  similar.  We  see  that  average  changes  in  the  wind  field  are  2  m  s'3-  in  the 
upper  troposphere  and  1.0  to  1.5  m  s"1  in  the  middle  and  lower  troposphere. 
The  non-divergent  wind  and  analyzed  temperatures  on  the  14  standard  pressure 
levels  are  then  interpolated  in  the  vertical  to  the  ten  model  sigma  levels. 
Our  diagnostic  equation  for  the  geopotential  is  then  used  to  compute  a 
balanced  geopotential  on  the  sigma  surfaces.  In  Fig.  3b  we  show  the  rms 
changes  in  the  temperature  at  each  of  the  model  sigma  levels  resulting  from 
solving  for  the  non-linear  mass  balance.  Changes  in  the  temperature  in  the 
upper  troposphere  of  1°C  reduce  to  0.5°C  in  the  middle  troposphere.  However  in 
the  lower  troposphere  large  changes  of  4.0°C  are  found.  For  comparison,  the 
actual  analyzed  winds  and  temperatures  are  interpolated  to  the  model  sigma 
levels.  In  Figure  4,  we  compare  the  wind  and  temperature  fields  for  the 
analyzed  fields  and  the  statically  initialized  fields,  after  having 
interpolated  them  back  to  500  and  1000  mb  pressure  levels  for  display.  At  the 
500  mb  level,  the  temperature  changes  are  small  although  the  small  trough  is 
weaker  as  a  result  of  the  non-linear  mass  balance.  At  the  1000  mb  level 
however  the  strength  of  the  front  is  greatly  weakened  by  the  scheme,  as  is 
reflected  in  the  large  rms  change  in  temperature  seen  in  Fig.  3b. 


4.  IMPLICIT  NORMAL  MODE  INITIALIZATION 


We  follow  the  vertical  mode  initialization  scheme  of  Bourke  and  McGregor 
(1983)  by  using  the  model  dynamical  equations  to  filter  out  the  high  frequency 
inertia-  gravity  waves.  Filtering  conditions  are  applied  to  the  model 
dynamical  equations  to  derive  linear  diagnostic  equations  for  the  mass 
divergence  and  generalized  geopotential,  which  are  solved  iteratively  for  the 
first  three  vertical  modes  of  the  numerical  model.  The  further  condition  that 
the  linearized  potential  vorticity  is  unchanged  by  the  procedure  is  required 
to  compute  the  vorticity.  Boundary  conditions  on  the  generalized  geopotential 
and  mass  divergence  are  required.  As  is  customary,  we  keep  the  geopotential, 
temperature  and  pressure  fixed  at  the  lateral  boundaries.  To  provide  a 
boundary  condition  for  the  divergence,  an  approximate  divergence  is  computed 
using  the  thermodynamic  equation.  In  our  scheme,  changes  in  the  tangential 
wind  along  the  lateral  boundaries  are  consistent  with  the  changes  in  the 
vorticity  and  divergence  computed  by  the  scheme. 

4.1  Inertia-  Gravity  Wave  Modes 


We  express  the  dynamical  equations  in  terms  of  the  time  tendencies  of  mass 
weighted  vorticity,  divergence  and  generalized  geopotential,  linearized  about 
the  basic  state  at  rest  with  mean  temperature  T*.  The  equations,  with  the  p 
term  included  with  the  non-linear  terms  on  the  RHS,  are 
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Here  f  is  the  Coriolis  parameter,  which  is  a  function  of  latitude  and  we  have 
further  defined  the  mass  weighted  vorticity  £  on  the  sigma  surfaces  as 
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In  the  equations  we  have  ignored  the  staggering  of  the  vorticity  and  mass 
divergence  relative  to  the  generalized  geopotential. 


The  freely  propagating  inertia-  gravity  waves  in  the  model  are  solutions 
to  the  linearized  equations  (4.1),  (4.2)  and  (4.3)  with  the  forcing  terms  on 
the  RHS  of  the  equations  equal  to  zero.  Taking  the  time  derivative  of  Eq. 
(4.2)  and  eliminating  the  tendencies  of  the  generalized  geopotential  and  the 
vorticity  using  Eqs.  (4.3)  and  (4.1)  with  the  terms  on  the  RHS  zero,  we  find 
that  the  mass  divergence  satisfies  the  equation 

92  D 
8t2 


f2  - 


(4.5) 


In  this  case,  the  vertical  modes,  which  are  again  found  by  solving  for  the 
eigenmodes  of  matrix  M3,  have  the  same  vertical  structure  as  that  computed  for 
the  gravity  wave  modes  in  section  2.1.  The  "ageostrophic  deviations" 
f£  -  V2*.  which  is  f  times  the  ageostrophic  vorticity,  can  also  be  shown  to 
satisfy  the  same  equation.  Eliminating  D  from  the  homogeneous  form  of  Eqs. 
(4.1)  and  (4.3),  an  expression  for  the  tendency  of  the  linearized  potential 
vorticity  Q  is  obtained 

§£  (  M3£  -  f  ♦  )  -  0  (4.6) 

where  the  linearized  potential  vorticity  Q  is  defined  as 

Q  =  (  £  -  f  M"1  ♦  )  (4.7) 


That  is,  the  linearized  potential  vorticity  Q  is  unchanged  by  the  inertia- 
gravity  wave  motions  (see  also  Errico,  1986). 

Now  in  terms  of  the  amplitude  of  the  mass  divergence  for  each  vertical 
mode,  we  have 
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By  including  the  Coriolis  parameter  in  the  linearization,  we  have  introduced  a 
low  frequency  cutoff  for  the  inertia-  gravity  wave  modes  (see  also  Gill,  1982, 
for  example).  Assuming  that  the  latitudinal  wavelength  of  the  modes  are  small 
compared  to  the  change  in  the  Coriolis  parameter  f  with  latitude,  Eq.  (4.8)  is 
a  wave  equation  with  the  phase  speed  of  the  plane  wave  modes  given  by 
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where  R<j  is  the  Rossby  radius  of  deformation  defined  by 


Rd 

(4.10) 

f 

and  ft  = 

(*x,  Ky)  is  the  horizontal  wavenumber  with  K2 

=  Kx2  +  Ky2.  At  the 

higher  frequencies  (that  is,  for  wavelengths  small  compared  to  2  t  R<j  ) ,  the 
phase  speeds  of  the  modes  can  be  approximately  given  by  and  the  modes 
are  essentially  the  gravity  wave  modes  computed  in  section  2.1.  For  our  ten 
xayer  model  and  our  limited  area  domains,  the  gravity  waves  for  the  first  two 
modes  are  largely  of  high  frequencies,  since  the  Rossby  radius  of  deformation 
is  large  (7,256  km  and  2,907  km  at  35°N  for  example,  for  the  first  and  second 
modes,  respectively),  compared  to  the  size  of  the  model  domain.  However,  for 
the  higher  modes,  the  Rossby  radius  of  deformation  becomes  comparable  or 
smaller  than  the  domain  size  (1,126  km  and  693  km  for  the  third  and  fourth 
modes  at  35°N,  for  example),  and  the  gravity  waves  are  of  lower  frequency. 

4 . 2  Vertical  Mode  Initialization 


4.2.1  Filtering  Equations.  We  want  to  reduce  the  initial  amplitude  and 
tendency  of  the  high  frequency  inertia-  gravity  waves  so  that  their  amplitude 
remains  small  during  the  integration  of  the  numerical  model.  By  requiring  that 
initially  the  first  and  second  derivatives  of  the  mass  divergence  are  zero, 
the  amplitudes  of  the  high  frequency  inertia-  gravity  waves  will  be  kept  small 
during  integration  of  the  model.  Taking  the  time  derivative  of  Eq.  (4.2)  and 
substituting  for  the  tendencies  of  the  generalized  geopotential  and  vorticity 
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from  Eqs.  (4.1)  and  (4.3),  an  expression  for  the  second  derivative  of  the  mass 
divergence  with  respect  to  time  is  obtained, 


82D 
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v2  a4  + 


(4.11) 


The  terms  on  RHS  of  Eqs.  (4.1),  (4.2),  (4.3),  which  include  the  beta  and  non¬ 
linear  terms,  vary  slowly  with  time  compared  to  the  time  scale  for  high 
frequency  gravity  waves.  Therefore  we  can  ignore  the  time  tendency  of  Ad  in 
Eq.  (4.11).  Then  applying  our  filtering  conditions 


8d  =  a^D 
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(4.12) 


to  Eqs.  (4.2)  and  (4.11)  respectively,  we  obtain 

V2  *  -  f  $■  -  Ajj  (4.13) 

M3  V2  D  -  f2  D  *  V2A}-  fAf  (4.14) 

To  complete  the  set  of  equations  a  further  condition  is  required.  Assuming 

that  the  changes  to  our  initial  fields  due  to  our  filtering  procedure  will  be 
small,  then  the  changes  to  our  fields  represent  that  part  of  the  fields  due  to 
the  freely  propagating  inertia-  gravity  waves.  If  they  are  small,  their  motion 
can  be  described  by  the  linearized  equations  (4.1),  (4.2)  and  (4.3)  with  the 
terms  on  the  RHS  equal  to  zero.  We  can  therefore  require  that  the  linearized 
potential  vorticity  be  unchanged  by  our  initialization  procedure.  Our 
filtering  conditions  are  now 

^  V2  ♦  -  f2  ♦  -  ^(Ap+fQjj)  (4.15) 

V2  D  -  f2  D  -  V2At-  fA^  (4.16) 

£  -  f  •  ■  Qq  ,  a  constant  (4.17) 


We  only  need  to  filter  the  high  frequency  gravity  waves,  whose  phase  speeds 
are  much  larger  than  the  typical  speeds  of  weather  systems.  Therefore  our 
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filtering  conditions  need  only  be  applied  for  those  vertical  modes  whose  phase 
speeds  are  larger  than  about  25  m  s"1.  For  the  ten  layer  model  example  shown 
in  Table  2,  filtering  the  first  three  modes  is  sufficient.  Using  the  same 
vertical  structure  we  defined  earlier  we  can  express  our  variables  in  terms  of 
amplitudes  for  each  of  the  vertical  modes,  where 
e  -  E-1#,  d  -  E_1D,  V  =  E"1^,  1)  «=  E-^-Q. 

Then  our  filtering  conditions  for  the  first  three  modes  k  =  1,  2,  3  can  be 
written 


V2  e 


k 


V, 


L. 

X,  k 
k 


(A. 18) 

(A. 19) 
(A. 20) 


where  the  forcing  functions  a^  and  b^  are  the  kfc^  elements  of  the  vectors 


E  <  +  f  %T7) 

e-1  (  v2  a,  . 


(A. 21) 

(A. 22) 


Here  we  have  reintroduced  the  averaging  operators  to  take  into  account  the 
staggering  of  the  variables  on  our  model  grid.  Since  the  forcing  functions  on 
the  RHS  of  Eqs.  (4.18)  and  (4.19)  depend  on  the  vorticity,  divergence  and 
geopotential,  the  set  of  equations  is  solved  iteratively.  With  a  first 
estimate  for  the  variables  given  by  the  uninitialized  fields,  the  forcing 
functions  for  the  Helmholtz  equations  can  be  computed  by  integrating  the  model 
to  obtain  the  non-linear  terms  Ap,  A$,  Ag,  and  the  initial  potential  vorticity 
Qq  computed.  By  solving  the  Helmholtz  equations  for  the  amplitudes  of 
generalized  geopotential  and  divergence,  new  values  of  the  variables  can  be 
computed.  By  recomputing  the  forcing  functions  the  process  can  be  repeated. 


We  require  boundary  conditions  for  the  amplitudes  of  the  generalized 
geopotential  and  the  mass  divergence  to  solve  the  Helmholtz  equations  at  each 
iteration  of  our  initialization  procedure.  We  choose  to  keep  ♦  fixed  at  the 
boundary,  so  that  at  the  boundary  T  and  ps  are  unchanged  by  the 
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Initialization  procedure.  It  has  been  customary  to  set  the  amplitude  of  the 

mass  divergence  in  the  first  three  modes  to  zero  along  the  boundary  and 

adjusting  the  integrated  divergence  over  the  model  domain.  However  for  a 

domain  with  substantial  topographic  features  along  the  boundary,  this  is  too 

restrictive.  We  choose  to  estimate  the  divergence  at  the  first  mass  point  in 

from  the  boundary  using  the  thermodynamic  equation  and  neglecting  the  tendency 

of  the  generalized  geopotential.  In  terms  of  the  amplitudes  of  the  vertical 

9e 

modes  we  have,  using  thermodynamic  equation  (4.3)  with  =  0, 
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4.2.2  Iterative  Procedure.  By  defining  incremental  changes  to  the 
amplitudes  of  mass  divergence,  generalized  geopotential  and  vorticity  for  each 
iteration  i 

"k  -  d£‘l  +  *4  •  4  -  'i'1  *  -  “k'1  +  K 


an  equivalent  scheme  can  be  derived  for  the  incremental  changes.  For  the  it^ 
iteration  of  the  amplitudes  of  the  generalized  geopotential,  mass  divergence 
and  vorticity,  we  have  for  the  ktlx  mode 
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(4.26) 


(4.27) 


We  can  use  Eqs.  (4.2)  and  (4.11)  to  compute  the  residuals  which  remained  after 
the  previous  (i-l)th  iteration.  Using  Eq.  (4.17)  to  substitute  for  the 
vorticity  in  the  equation  for  the  mass  divergence  Eq.  (4.2)  and  multiplying  by 
E"l,  we  find  that 
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That  is,  the  residual  remaining  after  the  (i-l)^  iteration  for  the  equation 
of  the  amplitude  of  the  generalized  geopotential  is  given  by  the  tendency  of 
the  mass  divergence.  A  similar  residual  can  be  computed  for  equation  for  the 
amplitude  of  the  mass  divergence.  Taking  the  time  derivative  of  Eq.  (4.2)  and 
substituting  for  the  second  time  derivative  of  the  mass  divergence  in  Eq. 
(4.11)  and  then  multiplying  by  E"1,  we  find 
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That  is,  the  residual  remaining  after  the  ( i- 1 ) iteration  is  given  by  the 
tendency  of  the  amplitude  of  the  ageostrophic  deviation  f£  -  V2#.  Then  by 
subtraction,  we  find  that  the  incremental  changes  are  forced  by 


(4.30) 
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where  the  terms  on  the  RHS  of  Eqs.  (4.30),  (4.31)  are  now  the  residuals  from 
the  previous  iteration.  To  compute  the  residuals  on  the  RHS,  we  integrate  the 
model  one  time  step  to  compute  adiabatic  tendencies  without  friction,  diabatic 
heating,  or  updating  of  the  values  of  the  model  variables  at  the  lateral 
boundaries.  The  residuals  are  then  computed  from 


0d 

at 


dv 

at 


8e 

at 


-1  0D 
E  8t 


E 


■  l  8$ 
8t 


-!  8* 
E  8t 


a^u 

r  b^v  i 

• 

e"1 

6 

*s 

+  £  6 

h7 

- 

X 

.  8t  . 

hx  7 

x  8t 

. 

E'1 

6 

*s 

1 

-  —  6 

Bp*  u 
h 

* 

X 

,  Bt  . 

h7  7 

x 

l  x  at 

. 

E'1 

“i 

'  BPJ  ■ 

+  [  R  T* 

-#*] 

.  8t  . 

1  8t 

J 

(4.33) 

(4.34) 

(4.35) 


27 


In  terms  of  the  changes  to  the  amplitude  of  the  generalized  geopotential  our 
boundary  condition  on  changes  to  e  becomes 

Aej^  *  0.  (4.36) 
In  terms  of  changes  to  the  amplitude  of  the  divergence  the  boundary  condition 
becomes 


(4.37) 


4.2.3  Changes  to  the  Horizontal  Wind.  Temperature  and  Surface  Pressure. 
The  resulting  changes  to  the  horizontal  wind  field  can  then  be  computed  from 
the  changes  in  the  divergence  and  vorticity  by  solving  for  the  changes  in  the 
velocity  potential  and  the  stream  function.  The  grid  stencil  for  the 
divergence  (and  velocity  potential)  and  the  vorticity  (and  stream  function)  is 
shown  in  Fig.  5.  That  is,  given  the  change  in  the  divergence  AD  *  E  Ad  we 
solve  a  Poisson’s  equation 

V2  A*  *  AD  (4.38) 
for  the  change  AX  in  the  velocity  potential  j  at  the  interior  mass 
(generalized  geopotential  •)  points  (see  Fig.  5)  for  each  sigma  level  of  the 
model.  The  boundary  condition  on  AX  is 

Ajr  -  0  (4.39) 
at  the  model  lateral  boundaries.  This  produces  no  change  in  the  tangential 
wind  along  the  model  boundaries,  but  only  change  to  the  divergent  wind  in  the 
interior.  Now  given  the  change  in  the  vorticity  Af  ■  E  Av  we  solve  a 
Poisson’s  equation 

V2  A*  -  Af  (4.40) 
for  the  change  At  in  the  stream  function  t  for  each  sigma  level  of  the  model. 
The  vorticity  changes  calculated  from  Eq.  (4.32)  are  specified  at  interior 
points  staggered  half  a  grid  distance  from  the  mass  and  wind  points.  Since  the 
integrated  vorticity  over  the  domain  may  change,  the  integrated  tangential 
wind  along  our  model  boundary  may  change.  We  therefore  prescribe  a  boundary 
condition  of  no  change  in  the  stream  function 

A#  «  0  (4.41) 
at  fictitious  boundary  staggered  half  a  grid  point  outside  our  model  boundary, 
implying  that  the  non-divergent  wind  does  not  change  there.  However  the  non- 
divergent  wind  and  hence  the  tangential  wind  does  change  along  the  model 
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boundary*  The  gradients  of  the  changes  of  the  stream  function  and  the  velocity 
potential  at  each  sigma  level  of  the  model  then  define  the  changes 
to  the  mass  weighted  wind, 

A(pfu)  =  6X  A}f  -  6y  At  (4.42) 

A(PsV)  ■  Sy  Lx  +  Sx  At  (4.43) 

With  the  changes  to  the  wind  field  defined  in  this  way,  the  changes  to  the 

wind  field  along  the  lateral  boundaries  of  the  model  are  consistent  with  the 

changes  in  the  vorticity  and  mass  divergence  over  the  model  domain. 


Following  Temperton  (1984),  the  changes  in  the  surface  pressure  p8,  and 
the  temperature  T  can  be  derived  directly  from  changes  in  the  generalized 
geopotential  4.  We  consider  the  linearized  equations  for  the  motion  of  the 
freely  propagating  gravity  waves 
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Now  eliminating  D  we  can  relate  the  tendencies  of  surface  pressure  and 
temperature  to  that  of  the  generalized  geopotential  for  gravity  wave  motions, 
so  that 
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Since  the  changes  derived  from  our  initialization  procedure  represent  the 
gravity  wave  part  of  the  flow,  we  can  assume  a  wave  solution  for  the  changes. 
The  changes  to  the  surface  pressure  and  temperature  are  then  related  to  the 
changes  in  the  generalized  geopotential  by 

Apg  -  N2T  M^1  A4  -  N2T  M^1  E  Ae  (4.48) 

A(pgT)  -  Mj  M^1  A4  -  M^1  E  Ae  (4.49) 
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4.3  Convergence  of  the  Vertical  Mode  Scheme 

To  test  our  vertical  mode  initialization  procedure,  we  use  the  12  hourly 
NMC  2.5°  hemispheric  analyses  for  the  period  12Z  January  23  to  12Z  January  29, 
1986,  for  the  period  of  the  second  Intensive  Observing  Period  (IOP)  of  GALE. 
The  initial  synoptic  situation,  showing  a  cold  front  moving  off  the  east  coast 
of  the  U.S.  and  a  low  in  the  Gulf  of  Alaska,  is  shown  in  Figs.  4  and  8.  During 
this  IOP,  a  coastal  front  develops  along  the  east  coast  of  the  U.S.  and 
subsequently,  a  cyclone  develops  offshore  when  the  low  system  from  the  Gulf 
of  Alaska  reaches  the  east  coast.  The  thirteen  analyses  for  the  period  are 
interpolated  to  the  model  coordinates  for  the  two  different  horizontal  grids 
with  differing  domain  size  and  resolution  and  then  initialized  with  the 
vertical  mode  scheme  for  the  first  three  vertical  modes  only.  The  US  grid 
covers  a  domain  including  the  continental  U.S. A.  with  a  horizontal  resolution 
of  2°  longitude  by  1.5°  latitude.  The  GALE  grid  covers  a  smaller  domain 
including  the  eastern  half  of  the  U.S. A.  and  extending  out  over  the  Atlantic 
to  52 . 5®W  with  a  finer  horizontal  resolution  of  0.5®  in  longitude  and 
latitude.  In  the  vertical,  both  grids  use  ten  equally  spaced  sigma  levels.  The 
smoothed  model  topography  used  for  each  grid  is  shown  in  Fig.  6a  and  b.  For 
the  GALE  grid,  a  case  with  smoother  topography  along  the  lateral  boundary, 
shown  in  Fig.  6c,  is  also  generated  by  merging  the  courser  topography  from  the 
US  grid  (Fig.  6a)  with  the  GALE  topography  (Fig.  6b)  in  a  boundary  zone  with  a 
width  of  five  degrees  in  latitude  and  longitude.  The  merging  is  carried  out  by 
linearly  interpolating  the  courser  topography  from  the  US  grid  to  the  GALE 
grid,  and  replacing  the  GALE  topography  at  the  first  seven  points  in  from  the 
boundary.  At  the  eight  to  tenth  points,  a  linear  combination  is  used  with 
weights  given  by  (0.75,  0.25),  (0.5,  0.5)  and  (0.25,  0.75)  for  the  course 
topography  and  GALE  topography,  respectively. 


4.3.1  Using  the  low  resolution  US  grid.  At  the  analysis  times,  the 
initial  amplitudes  for  the  mass  divergence  and  vorticity  are  computed  on  the 
US  grid  for  each  of  the  vertical  modes,  prior  to  initialization.  The  mean 
amplitudes  of  the  mass  divergence  and  vorticity  for  each  of  the  vertical  modes 
are  averaged  over  the  US  domain  and  in  time  over  the  period  of  interest  for 
the  thirteen  analyses  and  shown  in  Table  3.  The  mean  amplitude  of  the 
vorticity  decreases  with  increasing  mode  number,  while  the  mean  amplitude  of 
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TABLE  3:  The  mean  mass  divergence  and  vorticity  on  the 
sigma  surfaces  of  the  model  for  the  US  grid, 
averaged  for  the  week  of  12Z  January  23  to  122 
January  29,  1986. 


Mode  No. 

Mass  Divergence 
(dynes  cm"^  s"^-) 

Vorticity 
(dynes  cm“2  s"^-) 

1 

7.68 

82.04 

2 

7.54 

36.28 

3 

12.63 

36.98 

4 

9.08 

32.86 

5 

7.63 

17.53 

6 

6.32 

11.82 

7 

4.81 

8.14 

6 

3.81 

5.56 

9 

2.92 

4.41 

10 

1.82 

2.65 

the  mass  divergence  maximizes  at  the  third  mode.  As  expected,  the  amplitudes 
of  the  divergence  are  an  order  of  magnitude  smaller  than  the  amplitude  of 
vorticity  for  the  external  mode.  The  analyses,  interpolated  to  the  model  sigma 
coordinates  for  the  US  grid,  are  initialized  with  the  vertical  mode  scheme  for 
the  first  three  vertical  modes  only.  After  each  iteration  of  our  vertical  mode 
initialization  procedure,  the  resulting  root-mean- square  (rms)  changes  in  the 
amplitudes  of  the  mass  divergence,  vorticity,  generalized  geopotential  and 
surface  pressure  were  computed.  In  Fig.  7  we  show  these  rms  changes,  averaged 
for  all  the  analysis  times,  at  each  iteration  of  the  initialization  procedure. 
For  each  of  the  first  three  vertical  modes  initialized,  it  can  be  noted  that 
the  mean  rms  changes  in  the  amplitudes  of  the  divergence  at  the  first 
iteration  of  the  procedure  are  as  large  as  the  initial  amplitudes  themselves, 
while  the  changes  in  the  amplitude  of  the  vorticity  are  very  small  compared  to 
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their  initial  values.  For  the  first  two  modes,  the  changes  in  the  amplitude  of 
the  mass  divergence  decrease  rapidly  with  increasing  iterations,  with  the 
changes  being  very  small  after  just  two  iterations.  For  the  third  mode,  the 
changes  in  amplitudes  of  the  mass  divergence  do  not  decrease  as  rapidly.  The 
mean  rms  changes  in  the  surface  pressure,  shown  in  Fig.  7d  are  less  than  a  mb 
for  each  iteration. 

To  demonstrate  the  effect  of  the  number  of  modes  initialized  on  the 
changes  in  the  amplitude  of  the  mass  divergence  with  each  iteration,  the 
number  of  modes  initialized  is  varied  from  one  to  six  modes  for  the  case  of 
one  analysis  on  12Z  January  23,  1986.  In  Fig.  7e  we  show  the  changes  in  the 
mass  divergence  for  the  highest  mode  number  initialized  for  each  of  the  cases. 
The  rate  of  decrease  in  the  mass  divergence  changes  increases  as  the  higher 
order  modes  are  initialized,  and  in  fact  increases  with  iteration  for  the 
sixth  mode  initialized.  For  the  first  two  modes,  the  gravity  waves  are 

essentially  of  high  frequency,  since  the  Rossby  radius  of  deformation  is  so 
large  for  these  modes  (see  section  4.1).  As  the  mode  number  increases  however, 
the  Rossby  radius  of  deformation  decreases  and  the  frequency  of  the  inertia- 
gravity  v..:de  decreases  for  the  same  wavelengths.  The  tendencies  of  the  mass 
divergence  for  the  gravity  modes  are  then  much  less  for  the  higher  order 
modes,  and  convergence  of  the  scheme  would  be  expected  to  be  slower. 

As  an  illustration,  we  demonstrate  the  result  of  the  vertical  mode 

initialization  using  the  NMC  analysis  for  122  January  23,  1986.  The  initial 
analyzed  sea-level  pressure  and  wind  field  at  the  sigma  level  a  *  0.25  are 
shown  in  Figs.  8a  and  8b.  A  deep  surface  low  of  987  mb  lies  west  of  Greenland 
and  high  pressure  dominates  the  eastern  U.S.  Aloft  at  the  jet  level,  strong 
jet  maxima  of  about  55  and  48  m  s-1  straddle  a  trough  over  the  eastern  U.S., 
with  a  further  jet  maximum  upstream,  entering  the  domain  from  the  Gulf  of 
Alaska.  The  surface  pressure  change  resulting  from  the  vertical  mode 
initialization  is  shown  in  Fig.  8c.  The  surface  pressure  adjustments  due  to 
initialization  are  small,  with  an  rms  change  in  this  case  of  0.9  mb  and  at 
most  several  mb  in  places.  After  the  initialization,  the  vertical  motion  o 
in  the  middle  troposphere  at  a  sigma  level  0  *  0.45,  is  shown  in  Figure  8d. 

Upward  motion  relative  to  the  sigma  surface  (negative  values  of  the  vertical 

motion)  is  found  on  the  south  side  of  the  jet  off  the  coast  of  North  America 
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and  also  ahead  of  the  exit  region  of  the  jet  streak.  Sinking  motion  is  found 
in  the  region  of  the  high  surface  pressure  over  the  eastern  D.S.  In  sigma 
coordinates,  the  sigma  surfaces  follow  the  sloping  topography,  with  the  result 
that  westerly  flow  in  the  lee  of  mountains  gives  rise  to  strong  upward  motions 
on  the  sigma  surface.  Such  strong  topographic  signatures  are  seen  in  the  lee 
of  the  Rockies  and  the  Appalachian  mountains.  In  Fig.  9  we  compare  the 
contributions  to  the  vertical  motion  from  the  first  three  modes  (Fig.  9a), 
which  are  initialized,  and  from  the  remaining  modes  (Fig.  9b),  which  are 
unchanged  by  the  initialization  procedure.  The  strong  signal  due  to  the 
mountains  clearly  dominates  the  vertical  motion  computed  for  the  first  three 
modes,  while  smaller  synoptic  scale  magnitudes  are  apparent  in  the  remaining 
modes . 

The  vertical  mode  initialization  scheme  removes  that  portion  of  the 
initial  wind  and  mass  fields  that  describe  the  inertia-  gravity  waves  (for 
modes  1  to  3).  Such  structures  then  should  be  seen  in  the  changes  made  to  the 
analyzed  fields  resulting  from  the  initialization.  In  Fig.  10a  we  show  the 
wind  and  geopotential  height  changes  at  the  jet  level  a  *  0.25.  At  the  jet 
level,  we  see  that  the  geopotential  changes  reach  30  to  AO  gpm  in  places.  The 
wind  field  shows  flow  of  several  m  s”1  crossing  the  contours  of  geopotential 
height,  indicative  of  inertia-  gravity  wave  structures  (see  Matsuno,  1966,  for 
example).  At  the  same  level  the  changes  in  temperature  and  the  u  and  v 
components  of  the  wind  are  shown  in  Fig.  10b,  c  and  d.  The  changes  in  the 
temperature  are  at  most  a  degree,  while  changes  in  the  wind  components  are  at 
most  several  m  6-1.  Typical  vertical  profiles  of  the  rms  changes  in  the  u 
component  and  the  temperature  are  shown  in  Fig.  11.  The  resulting  changes  in 
the  mean  temperature,  shown  in  Table  A,  are  very  small  and  at  most  0.15#C  in 
the  upper  troposphere. 

A. 3. 2  Using  the  high  resolution  GALE  grid.  For  the  smaller  GALE  grid  with 
0.5  degree  resolution,  the  mean  amplitudes  for  the  vorticity  and  mass 
divergence  are  shown  in  Table  5,  for  each  of  the  vertical  modes.  For  each 
iteration  of  our  vertical  mode  initialization  scheme,  the  mean  rms  changes  in 
the  amplitudes  of  the  divergence  and  the  vorticity  are  shown  in  Fig.  12,  for 
the  first  three  modes  initialized.  In  this  case,  the  amplitude  changes  for  the 
geopotential  (or  equivalently  vorticity)  decrease  very  rapidly  for  all  the 
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TABLE  4:  The  root-mean- Gquare  changes  in  the  mean 

temperature  on  the  US  grid  for  12Z  January  23, 
1986,  resulting  from  vertical  mode  initialization. 


Model  Level 

Sigma  Level 

AT* 

(°K) 

1 

0.05 

0.052 

2 

0.15 

0.143 

3 

0.25 

0.098 

4 

0.35 

0.048 

5 

0.45 

0.022 

6 

0.55 

0.018 

7 

0.65 

0.014 

8 

0.75 

0.069 

9 

0.85 

-0.001 

10 

0.95 

-0.011 

three  modes.  However,  the  amplitude  changes  for  the  mass  divergence  in  the 
case  of  the  third  mode  do  not  decrease  as  rapidly  and  in  fact  do  not  reduce  to 
zero.  On  the  GALE  grid,  smaller  scale  topographic  features  of  appreciable 
amplitude  are  present  (see  Fig.  6b),  compared  to  the  US  grid  (see  Fig.  6a). 
For  the  third  mode,  the  gravity  modes  at  these  smaller  scales  are  of  high 
frequency  and  of  shorter  vertical  scale.  Non-linear  effects  of  the  gravity 
waves  interacting  with  the  topography  can  become  significant.  Also  in  the 
boundary  zone,  where  there  was  no  special  smoothing  of  the  topography,  the 
divergence  forced  by  the  topography  for  these  modes  can  get  quite  large.  Then 
the  errors  in  the  computed  boundary  divergence  can  be  large  also.  To  test  the 
latter,  the  initialization  is  repeated  for  two  cases  with  the  topography  on 
the  GALE  grid  smoothed  in  a  five  degree  boundary  zone,  as  outlined  in  the 
beginning  of  section  4d  (see  Fig.  6c).  In  Fig.  12c  we  compare  the  changes  in 
the  amplitude  of  the  mass  divergence  when  the  unsmoothed  (Fig.  6b)  and  the 
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smoothed  topography  (Fig.  6c)  are  used  for  a  typical  example  (such  as  for  the 

analysis  on  12Z  January  23,  1986)  and  for  an  extreme  example  with  much  larger 

changes  at  higher  iterations  (for  122  January  24).  In  both  cases,  the  changes 
in  the  mass  divergence  decrease  much  more  rapidly  when  the  topography  is 
smoothed  in  a  five  degree  boundary  zone. 

For  the  smaller  GALE  grid,  the  initial  sea  level  pressure  and  wind  field 
in  the  upper  troposphere  at  sigma  level  a  *  0.25  for  12Z  January  23,  1986  is 

shown  in  Fig.  13a  and  b.  A  front  is  shown  moving  off  the  east  coast  of  the 

U.S.  with  high  pressure  dominating  the  eastern  half  of  the  U.S.  A  strong  jet 
streak  is  leaving  the  domain  at  the  north-east  comer  of  the  domain  and  a 
minor  short  wave  trough  with  its  associated  jet  maximum  is  located  on  the  Gulf 
of  Mexico.  The  surface  pressure  changes,  shown  in  Fig.  13c  are  small,  being  at 

TABLE  5:  The  mean  mass  divergence  and  vorticity  on  the 
sigma  surfaces  of  the  model  for  the  GALE  grid, 
averaged  for  the  week  of  12Z  January  23  to  12Z 
January  29,  1986. 


Mode  No. 

Mass  Divergence 
(dynes  cm*2  s*1) 

Vorticity 
(dynes  cm*2  s"1) 

1 

7.53 

92.17 

2 

8.12 

39.21 

3 

15.11 

42.80 

4 

10.05 

37.65 

5 

9.17 

19.48 

6 

6.73 

12.05 

7 

4.79 

7.96 

8 

3.84 

5.72 

9 

2.74 

4.32 

10 

1.72 

2.51 

most  a  mb.  In  this  GALE  domain,  the  western  boundary  is  at  about  1  km.  sloping 
down  to  500  to  250  meters  within  5  degrees  in  from  the  boundary  (Fig.  6b).  The 
topographic  features  in  the  domain  produce  a  strong  signal  in  the  vertical 
motion  field  shown  in  Fig.  13d.  The  increasing  westerly  shear  with  height  at 
the  western  boundary  and  in  the  lee  of  the  Appalachians,  produces  strong 
rising  motion  in  sigma  coordinates.  Away  from  the  influences  of  sloping  sigma 
surfaces  over  the  topography,  sinking  motion  is  observed  in  the  high  pressure 
along  the  Mississippi  river  valley  and  rising  motion  on  the  south  side  of  the 
jet  axis  off  the  east  coast.  Qualitative  agreement  is  found  with  the  vertical 
motion  field  produced  for  the  larger  US  domain.  When  the  smoother  topography 
in  the  boundary  zone  (see  Fig.  6c)  is  used,  the  large  noisy  values  in  the 
vertical  motion  along  the  northern  boundary  and  south-west  corner  (seen  in 
Fig.  13d)  are  removed,  while  the  interior  remains  unchanged  (see  Fig.  13e). 
The  strong  signal  remaining  in  the  vertical  motion  field  along  the  western 
boundary,  due  to  the  high  gradients  in  the  topography  there,  demonstrates  the 
importance  of  the  specification  of  accurate  boundary  conditions  used  for  the 
solution  of  the  divergence  changes  in  the  initialization  procedure.  For  this 
case  the  resulting  changes  in  the  geopotential  and  wind  and  temperature 
changes  at  a  =0.25  are  shown  in  Fig  14.  In  the  upper  troposphere,  the 
resulting  rms  changes  in  the  wind  components  are  about  1  m  s-1,  while  the  rms 
temperature  changes  are  about  0.5°C.  In  the  vertical,  the  variation  of  the  rms 
changes  in  the  wind  components  and  temperature  with  each  iteration  is  similar 
to  that  shown  for  the  US  grid  in  Fig.  11,  decreasing  with  increasing  number  of 
iterations.  Three  iterations  are  sufficient  for  the  changes  in  the  wind  and 
temperature  to  be  very  small. 
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5.  MODEL  INTEGRATIONS  WITH  STATIC  INITIALIZATION 

To  test  the  effect  of  the  split-explicit  time  integration  scheme  in 
smoothing  the  unwanted  high  frequency  oscillations  in  the  NRL  model,  the  model 
is  integrated  on  the  GALE  grid  with  uninitialized  data  both  with  the  split- 
explicit  scheme  described  in  section  2.2  and  with  an  explicit  time  integration 
scheme  of  smaller  time  steps.  Integrations  with  varying  degrees  of  static 
initialization  are  then  compared  in  three  other  experiments.  Time  series  of 
the  surface  pressure  and  the  vertical  motion  0  in  sigma  coordinates  at 
selected  points  on  the  GALE  grid  were  compared  in  the  five  experiments,  which 
are  listed  in  Table  6.  The  NMC  2.5  degree  hemispheric  analysis  for  12Z  January 
23,  1986  is  used  to  start  the  integrations  for  each  of  the  experiments.  This 
is  a  case  of  a  cold-air  damming  and  coastal  front  event,  which  occurred  from 
January  23-25,  1986,  during  GALE.  The  GALE  grid  covers  a  domain  from  22.5°S  to 
47.5°N  in  latitude  and  102. 5°W  to  57.5°W  in  longitude,  with  0.5  degree 
horizontal  resolution. 

For  the  experiments  here,  the  Perkey  Kreitzberg  lateral  boundary 
formulation,  described  in  section  2.3a  is  used.  To  provide  the  boundary 
tendencies  for  each  of  the  experiments,  the  non-divergent  wind  is  computed  for 


TABLE  6s  Experiments  with  Static  Initialization  on  GALE  grid. 


1A  Explicit  Integration  with  uninitialized  initial  state  using 
leapfrog  scheme  with  a  2  At  time  step  of  60  secs. 

IB  Split-Explicit  Integration  with  uninitialized  initial  state, 
using  a  2  At  time  step  of  300  secs. 

1C  Split-Explicit  integration  with  Non-divergent  Initial  State, 
using  a  2  At  time  step  of  300  secs. 

ID  Split-Explicit  integration  with  Static  Non-linear  Mass  Balance, 
using  a  2  At  time  step  of  300  secs. 


Supplementary  experiment: 

1C*  Same  as  1C,  but  with  boundary  tendencies  computed  from  non- 
divergent  wind  and  observed  temperatures. 
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each  of  the  12  hourly  NMC  analyses,  the  fields  interpolated  to  the  sigma 
coordinates  of  the  model  and  a  statically  balanced  temperature  computed  as 
outlined  in  section  3.2.  Initial  boundary  values  and  12  hour  tendencies  are 
then  extracted  for  the  5  point  boundary  zone,  defined  in  section  2.3a.  To 
provide  the  same  boundary  values  in  each  of  the  experiments,  the  initial 
fields  for  each  experiment  are  merged  with  the  statically  balanced  fields  in 
the  boundary  zone,  using  the  linear  function  a  defined  by  Eq.  (2.19).  That  is, 
we  multiply  the  initial  fields  by  (1-a)  and  add  a  times  the  statically 
balanced  fields  in  the  boundary  zone. 

5.1  Damping  bv  Split-Explicit  Scheme 

In  expt  1A,  the  model  is  integrated  in  time  with  a  conventional  leapfrog 
scheme  for  12  hours  with  a  2  At  time  step  of  60  s  starting  from  uninitialized 
data.  Oscillations  of  surface  pressure  of  as  much  as  5  to  8  mb  of  amplitude 
and  periods  of  1  to  2  hours  are  observed  in  the  first  12  hours  of  integration. 
Curve  A  in  Fig.  15a  shows  these  typical  oscillations  in  the  surface  pressure 
at  a  grid  point  at  90°W  and  35°N  in  the  western  half  of  the  domain.  Curve  A  in 
Fig.  15b  shows  the  vertical  motion  0  interpolated  to  a  sigma  level  at  a  =  0.5, 
for  the  same  grid  point.  The  curve  A  shows  a  typical  rapid  adjustment 

(increase  in  this  figure)  in  the  first  6  of  the  integration  with  smaller 
oscillations  of  periods  of  2  to  4  hours  superimposed.  The  higher  frequency 
oscillations  in  surface  pressure  are  largely  due  to  the  barotropic  external 
gravity  mode  while  the  adjustment  and  oscillations  in  the  vertical  motion  in 
the  middle  troposphere  are  largely  due  to  the  internal  gravity  modes. 

In  expt  IB,  the  model  is  integrated  for  48  hours  with  the  split-explicit 
scheme  using  a  2  At  time  step  of  300  s  starting  again  from  the  same 

uninitialized  data.  For  the  first  hour  or  so  the  oscillations  in  the  surface 

pressure  (Curve  B  in  Fig.  15a)  are  the  same  as  in  the  explicit  integration. 
However  the  oscillations  are  strongly  damped  in  the  next  three  hours  of 
integration.  Little  difference  is  noticed  in  Fig.  15b  in  the  variation  of  a 
at  O  *0.5  between  the  explicit  (Curve  A)  and  split-explicit  (Curve  B) 
integrations.  One  can  conclude  that  the  split-explicit  integration  scheme  acts 
to  reduce  the  amplitude  of  the  unwanted  external  gravity  waves  in  the  first 
chrej  to  four  hours  of  integration. 
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5.2  Non-Divergent  and  Statically  Balanced  Initial  Fields 


As  outlined  in  section  3.1,  the  analyzed  vorticity  on  pressure  surfaces  is 
used  to  first  calculate  the  non-divergent  component  of  the  wind.  A  first  guess 
of  the  surface  pressure  and  temperature  is  found  by  interpolation  to  the  model 
topography.  The  non-divergent  wind,  analyzed  temperature  and  humidity  are  then 
interpolated  to  the  model  sigma  levels.  The  initial  and  non-divergent  wind 
fields  at  1000  mb  and  500  mb  were  compared  in  Fig.  4.  In  expt  1C,  a  48  hour 
integration  using  the  split-explicit  scheme  is  then  performed  with  this  data. 
As  shown  by  Curve  C  in  Figure  16a,  the  amplitude  of  the  initial  oscillations 
of  the  surface  pressure  are  reduced  to  2  to  3  mb  and  are  largely  damped  out 
after  3  hours.  Oscillations  in  the  surface  pressure  of  a  mb  can  be  still  seen 
at  a  grid  point  at  83°W  and  35°N  (not  shown)  over  the  Appalachian  mountains. 
Curve  C  in  Figure  16b  shows  that  the  strong  adjustment  (increase)  in  the 
vertical  motion  a  is  still  present  in  the  first  6  hours  and  the  superimposed 
higher  frequency  oscillations  are  only  slightly  reduced  in  amplitude.  By 
removing  the  divergent  component  of  the  wind  from  the  analyzed  data  the 
initial  value  of  the  vertical  motion  is  also  reduced.  On  sigma  surfaces  over 
sloping  topography,  a  vertical  shear  of  the  non-divergent  wind  will  introduce 
divergence  on  the  sigma  surfaces,  as  will  errors  caused  by  vertical 
interpolation.  By  largely  removing  the  horizontal  divergence  we  have 
essentially  removed  the  external  gravity  mode  after  three  hours  of 
integration. 

In  a  supplementary  experiment  1C’,  boundary  tendencies  are  computed  from 
analyzed  temperatures  and  the  non-divergent  winds  on  pressure  surfaces, 
instead  of  from  statically  balanced  temperatures.  For  initial  conditions,  the 
analyzed  temperatures  and  winds  are  interpolated  to  the  sigma  coordinates 
without  merging  the  statically  balanced  temperatures  in  the  boundary  zone.  The 
integrations  in  this  case  were  largely  indistinguishable  from  those  expt  1C. 
Since  in  the  PK  scheme,  we  damp  only  the  tendencies  in  the  boundary  zone, 
initial  differences  in  the  boundary  zone  in  the  two  cases  are  damped,  leading 
to  similar  integrations.  Little  difference  was  also  found  in  the  integrations 
if  the  mass  divergence  is  removed  from  the  initial  fields  on  the  sigma 
surfaces  instead  of  on  the  pressure  surfaces. 
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A  static  Initialization  of  the  mass  field  is  performed  using  the  non¬ 
linear  mas 8  balance  equation  in  sigma  coordinates,  as  outlined  in  section  3.2. 
The  non-linear  balance  equation  is  used  to  derive  a  balanced  geopotential  and 
temperature  field  from  the  non-divergent  wind  field.  In  Fig.  4  we  showed  the 
initial  and  balanced  temperature  fields  interpolated  to  the  1000  mb  pressure 
level.  As  mentioned  in  section  3.3,  the  balanced  temperature  field  is  much 
smoother  than  the  initial  field,  due  to  the  smoothing  of  the  Laplacian 
operator.  In  ezpt  ID,  a  48  hour  integration  was  carried  out  starting  from  this 
statically  initialized  data.  The  oscillations  in  the  surface  pressure,  shown 
in  Curve  D  of  Figure  16a,  are  of  the  same  amplitude  as  in  the  non-divergent 
case  (Curve  C  in  Fig.  16a),  being  damped  after  3  hours  of  integration. 
However,  a  small  mean  drift  of  about  a  mb  is  seen  to  develop  in  the  surface 
pressure  (curve  D)  during  the  12  hours  of  integration.  Again  oscillations  of 
surface  pressure  of  a  mb  still  remain  over  the  Appalachian  mountains.  The 
vertical  motion  shown  in  curve  D  of  Figure  16b  still  shows  the  rapid 
adjustment  (increase)  during  the  first  5  hours,  but  the  higher  frequency 
oscillations  are  much  reduced  in  amplitude  and  mostly  eliminated  after  4  hours 
of  integration.  The  non-linear  balance  of  the  mass  field  essentially  removes 
the  internal  gravity  waves  except  for  the  initial  adjustment  in  the  first  4 
hours . 

When  using  the  split-explicit  scheme  and  starting  from  uninitialized  or 
initialized  data,  the  12  to  48  hour  forecasts  are  very  similar.  Even  in  expt 
ID,  where  the  static  initialization  had  smoothed  out  the  initial  temperature 
gradient  along  the  front  (see  Fig.  4),  the  model  regenerates  the  temperature 
gradient  and  intensifies  it  further  in  the  first  12  hours.  As  an  example,  we 
show  in  Fig.  17,  the  variation  of  the  surface  pressure  and  the  vertical  motion 
at  a  grid  point  in  the  second  twelve  hours  of  integration  for  expts  IB,  1C  and 
ID.  In  the  case  of  expts  IB  and  1C,  the  variation  in  the  surface  pressure  is 
much  the  same.  Since  the  same  boundary  values  and  tendencies  were  used  in  the 
three  expts  (expts  IB,  1C  and  ID),  we  see  that  the  drift  in  the  surface 
pressure,  produced  by  the  static  initialization,  is  eliminated  in  the  second 
twelve  hours.  Only  in  the  case  of  expt  1A,  with  a  12  hour  explicit  time 
integration  with  uninitialized  data,  are  high  frequency  oscillations  in  the 
surface  pressure  and  vertical  motion  still  substantial  after  12  hours  of 
integration,  producing  more  noticeable  differences  in  the  12  hour  forecast. 
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6. 


MODEL  INTEGRATIONS  WITH  VERTICAL  MODE  INITIALIZATION 


To  assess  the  effect  of  the  vertical  mode  initialization  procedure  in 
removing  the  high  frequency  gravity  wave  oscillations,  the  NRL  model  is 
integrated  starting  from  initialized  and  uninitialized  data,  for  the  two 
different  US  and  GALE  grids  of  differing  domain  size  and  horizontal 
resolution.  The  influence  of  two  different  lateral  boundary  treatments,  namely 
the  tendency  relaxation  scheme  of  Perkey  and  Kreitzberg  (PK)  and  the  Davies 
relaxation  scheme,  is  also  investigated.  The  experiments  are  summarized  in 
Table  7.  As  described  in  section  2.3,  model  computed  tendencies  are  relaxed  to 
specified  boundary  tendencies  in  a  boundary  zone  of  5  points  at  each  time  step 
in  the  PK  scheme.  The  12  hourly  boundary  tendencies,  derived  from  the  NMC 
hemispheric  analyses  interpolated  to  the  model  grid.  In  the  Davies  scheme, 
large  scale  boundary  values  are  computed  by  linear  interpolation  from  12 
hourly  values,  derived  from  NMC  analyses,  which  have  been  interpolated  to  the 
model  grid.  The  computed  model  variables  are  then  relaxed  to  the  boundary 


TABLE  7:  Experiments  with  Vertical  Mode  Initialization. 

US  GRID  (with  a  2  At  time  step  of  400  secs): 

2B  24  Hour  integration  with  uninitialized  initial  state,  using 
PK  lateral  boundary  scheme. 

2C  24  Hour  integration  with  initialized  initial  state,  using 
PK  lateral  boundary  scheme. 

2E  24  Hour  integration  with  initialized  initial  state,  using 
Davies  lateral  boundary  scheme. 

GALE  GRID  (with  a  2  At  time  step  of  100  secs): 

3B  12  Hour  integration  with  uninitialized  initial  state,  using 
PK  lateral  boundary  scheme. 

3C  12  Hour  integration  with  initialized  initial  state,  using 
PK  lateral  boundary  scheme. 

3D  12  Hour  integration  with  uninitialized  initial  state,  using 
Davies  lateral  boundary  scheme. 

3E  12  Hour  integration  with  initialized  initial  state,  using 
Davies  lateral  boundary  scheme. 
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values  themselves  at  each  time  step  in  a  boundary  zone  of  6  points.  For 
integrations  with  initialized  fields,  the  analyzed  fields  interpolated  to  the 
model  grid  are  also  initialized  to  provide  balanced  boundary  values  and 
tendencies.  In  this  way  we  prevent  a  lot  of  noise  propagating  from  the 
boundary  into  the  interior  of  the  model  domain,  which  could  contaminate  the 
results.  On  the  two  grids,  the  smoothed  topography,  as  shown  in  Fig.  6  is  used 
for  the  integrations. 

6.1  Integrations  on  the  US  grid 

In  the  first  three  experiments  listed  in  Table  7,  we  compare  24  hour 
integrations  with  initialized  and  uninitialized  initial  conditions  on  the  US 
grid,  using  the  two  different  lateral  boundary  treatments.  The  split-explicit 
integration  scheme  is  used  in  each  case  with  a  2  At  time  step  of  400  seconds. 
In  expt  2B,  the  PK  tendency  relaxation  scheme  is  used  for  the  lateral  boundary 
treatment  with  uninitialized  initial  conditions.  The  12  hourly  boundary 
tendencies  are  derived  from  the  uninitialized  analyzed  fields.  Initial 
oscillations  of  nearly  8  mb  in  the  surface  pressure,  caused  by  the 
uninitialized  initial  conditions,  are  damped  by  the  split-explicit  scheme  in 
the  first  6  hours  or  so,  as  was  shown  in  section  5.1.  As  an  example,  curve  B 
in  Figs.  18a  and  b  shows  the  variation  of  the  surface  pressure  with  time  at  a 
grid  point  at  90°W  and  35°N  on  the  US  grid.  In  expt  2C,  initialized  initial 
conditions  are  used  with  the  PK  lateral  boundary  scheme.  The  boundary 
tendencies  in  this  case  are  derived  from  initialized  NMC  analysis  fields.  With 
initialized  initial  conditions,  the  high  frequency  oscillations  are  almost 
completely  removed,  even  over  high  topography  such  as  the  Rocky  Mountains.  An 
example  with  the  surface  pressure  can  be  seen  with  curve  C  in  Fig6.  18a  and  b. 
Large  high  frequency  oscillations  in  the  vertical  motion  are  also  eliminated 
for  the  interior  grid  points,  as  can  be  seen  in  Figs.  18c  and  d.  Even  though 
both  expts  2B  and  2C  used  somewhat  different  boundary  conditions 
(uninitialized  values  versus  initialized  values)  no  mean  drift  was  produced. 
In  the  boundary  zone,  the  difference  seen  between  the  results  in  expt  2B  and 
that  in  expt  2C,  is  due  to  gravity  waves,  which  are  damped  in  time  by  the  PK 
scheme. 
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In  expt  2E,  the  Davies  scheme  is  used  with  initialized  initial  conditions. 
As  can  be  seen  by  Curve  E  in  Fig.  18,  use  of  the  Davies  scheme  produces  a 
similarly  smooth  integration  of  the  surface  pressure  and  the  same  variation  in 
the  vertical  motion.  However,  the  variation  of  the  surface  pressure  with  time 
in  this  case  does  differ  somewhat  after  4  hours  from  that  found  with  the  PR 
tendency  relaxation  scheme.  The  PR  scheme  produced  some  artificially  larger 
vertical  motions  at  the  first  grid  point  inside  the  boundary.  The  time 
variation  of  the  surface  pressure  and  vertical  motion  at  a  grid  point  in  the 
boundary  zone,  at  a  distance  of  three  times  the  grid  spacing  in  from  the 
boundary,  is  shown  in  Fig.  19.  A  linear  variation  of  the  surface  pressure  is 
seen  with  time  in  the  boundary  zone  when  the  Davies  scheme  is  used.  With  the 
PR  scheme  oscillations  of  long  period  of  about  12  hours  can  be  seen  in  the 
surface  pressure  and  the  vertical  motion.  The  variation  of  the  vertical  motion 
is  small  for  the  Davies  scheme  at  this  point.  This  noise,  generated  by  the  PR 
scheme  in  the  boundary  zone,  may  have  propagated  into  the  interior  to  cause 
the  differences  seen  in  the  integrations  after  4  hours.  The  difference  can  be 
seen  (in  Fig.  18b  for  example)  as  a  low  amplitude  oscillation  of  long  period 
at  12  to  24  hours  integration.  However  little  overall  differences  can  be  seen 
in  the  12  and  24  hour  forecast  fields  between  each  of  the  experiments. 

6.2  intggmlang-aa  .the  gape  grid 

A  series  of  four  integrations  are  conducted  on  the  GALE  grid  with 
uninitialized  and  initialized  initial  conditions  with  the  two  lateral  boundary 
schemes  (see  last  four  experiments  listed  in  Table  7).  In  each  case,  the 
split-explicit  scheme  is  used  for  time  integration  with  a  2  At  time  step  of 
200  seconds.  For  the  cases  with  initialized  initial  conditions,  12  hourly 
boundary  values  and  tendencies  are  computed  from  the  NMC  analyzes, 
interpolated  to  the  model  grid  and  then  initialized.  The  integrations  with 
initialized  initial  conditions,  expts  3C  and  3E,  essentially  remove  the  high 
frequency  oscillations.  As  an  example  in  Fig.  20,  we  show  the  variation  of 
surface  pressure  and  the  vertical  motion  in  the  middle  troposphere  at  two  grid 
points  on  the  GALE  grid  for  the  different  integrations  with  initialized  and 
uninitialized  initial  conditions.  The  grid  point  at  90°W  and  35°N  is  in  the 
western  half  of  the  domain,  while  the  grid  point  at  83CW  and  35°N  lies  on  the 
top  of  the  Appalachian  Mountains.  With  the  PR  lateral  boundary  scheme. 
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oscillations  of  a  mb  or  so  are  seen  during  the  integration  period  when 
uninitialized  initial  conditions  are  used  (Curve  B).  With  initialized  initial 
conditions  obtained  with  the  vertical  mode  scheme,  the  oscillations  are 
removed  except  for  an  initial  hump  of  0.5  mb  in  the  surface  pressure  seen 
after  the  first  two  hours  of  integration  (Curve  C) .  With  the  Davies  scheme  and 
initialized  initial  conditions,  a  smoother  variation  of  the  surface  pressure 
is  seen  (Curve  E) .  The  rapid  initial  adjustment  in  the  vertical  motion  has 
been  largely  removed,  leaving  a  slower  increase  in  the  vertical  motion  with 
time. 

In  expt  3D,  the  Davies  scheme  for  the  lateral  boundary  treatment  is  used 
with  uninitialized  initial  conditions.  A  different  response  is  produced  in  the 
surface  pressure,  while  the  response  in  the  vertical  motion  field  is  similar 
to  the  PK  scheme  for  grid  points  in  the  interior  of  the  domain.  In  Fig.  21,  we 
compare  the  surface  pressure  variation  and  vertical  motion  in  the  middle 
troposphere  at  our  two  grid  points  for  initialized  and  uninitialized  initial 
conditions,  when  the  Davies  scheme  is  used.  In  the  case  of  the  Davies  scheme 
with  uninitialized  initial  conditions,  an  initial  shock  of  large  amplitude  is 
seen  in  the  surface  pressure  (curve  D  in  Figs.  21a  and  b),  which  is  rapidly 
damped  in  the  first  4  hours  or  so.  The  scheme  acts  to  damp  any  gravity  waves 
that  propagate  into  the  boundary  zone  from  the  interior.  Using  initialized 
initial  conditions,  this  initial  shock  is  eliminated  (curve  E  in  Figs.  21a  and 
b) .  With  time,  the  integrations  in  the  first  12  hours  differ  by  as  much  as  a 
mb  or  so  in  expts  3D  and  3E  (see  Fig.  21).  The  difference  is  explained  by  the 
fact  that  in  expt  3D,  we  force  uninitialized  boundary  values  in  the  boundary 
zone,  and  the  model  solution  in  the  interior  is  forced  to  adjust  in  the  first 
three  of  four  hours  of  integration. 

In  the  boundary  zone  large  differences  are  found  when  using  the  different 
boundary  treatments.  In  Figs.  22a  and  b,  we  show  the  variation  of  the  surface 
pressure  at  two  grid  points  in  the  boundary  zone  for  the  different  lateral 
boundary  treatments.  The  grid  point  at  59°W  and  40°N  lies  a  distance  of  three 
grid  lengths  from  the  lateral  boundary  (at  57.5®W),  while  the  grid  point  at 
58°W  and  40®N  is  a  distance  of  one  grid  length.  With  the  Davies  scheme,  the 
boundary  values  are  strongly  forced  in  the  boundary  zone,  which  is  reflected 
in  the  linear  variation  in  the  surface  pressure  in  the  boundary  zone  (see 
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curve  E  in  Figs.  22a  and  b,  for  example).  With  the  PK  scheme,  oscillations  in 
the  surface  pressure,  with  an  amplitude  of  one  to  two  mb  and  period  of  about 
twelve  hours,  can  be  seen  in  the  boundary  zone  (see  curve  C  in  Fig.  22a  and 
b) .  The  two  schemes  also  differ  in  their  response  in  the  vertical  motion  in 
the  boundary  zone.  With  the  PK  scheme,  the  vertical  motion  piles  up  at  the 
first  grid  point  in  from  the  boundary,  while  with  the  Davies  scheme,  larger 
vertical  motions  are  found  in  the  zone  between  the  interior  region  and  the 
boundary  zone.  The  variation  of  the  vertical  motion  at  the  two  grid  points  in 
the  boundary  zone  is  shown  in  Figs.  22c  and  d.  Using  the  Davies  scheme,  curve 
E  in  Fig.  22c  shows  a  slow  increase  of  the  vertical  motion  to  20  x  10"3  hr"*, 
for  the  grid  point  which  is  a  distance  of  three  grid  lengths  in  from  the 
lateral  boundary.  In  Fig.  22d,  curve  E  shows  a  slower  linear  variation  of  the 
vertical  motion,  increasing  to  10  x  10“3  hr"1,  at  the  grid  point  one  grid 
length  from  the  boundary.  With  the  PK  scheme  low  amplitude  changes  in  the 
vertical  motion  occur  at  the  boundary  between  the  boundary  zone  and  the 
interior  (see  curve  C  in  Fig.  22c),  while  unrealistically  large  values  of 
30  x  10“3  hr-1  are  reached  in  the  vertical  motion  at  the  grid  point  lying  a 
distance  of  one  grid  length  in  from  the  lateral  boundary  (see  curve  C  in  Fig. 
22d) .  In  Fig.  23,  we  compare  the  surface  pressure  variation  and  vertical 
motion  in  the  middle  troposphere  at  our  two  grid  points  in  the  boundary  zone 
for  initialized  and  uninitialized  initial  conditions,  when  the  Davies  scheme 
is  used.  The  boundary  values  are  again  strongly  forced  in  the  boundary  zone  in 
both  cases.  The  initial  shock  in  the  surface  pressure  variation  in  curve  D,  in 
Fig.  23a,  is  eliminated  by  the  initialization  (curve  E) . 
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7.  SUMMARY  AND  CONCLUSIONS 

To  remove  the  high  frequency  gravity  wave  oscillations,  various 
initialization  procedures  have  been  tested  for  use  with  the  Naval  Research 
Laboratory  limited  area  numerical  weather  prediction  model.  Operational 
analyses  obtained  from  the  National  Meteorological  Center  (NMC)  for  the  period 
of  the  second  intensive  observing  period  of  the  GALE  experiment  are  used  to 
test  the  procedures  and  provide  initial  and  boundary  conditions  for  model 
integrations.  The  model  is  integrated  with  initial  conditions  with  varying 
degrees  and  type  of  initialization  on  two  different  model  grids,  one  a  low 
resolution  grid  of  2°  longitude  by  1.5°  latitude  covering  the  continental  U.S. 
(US  grid)  and  the  other  a  higher  resolution  grid  of  0.5°  in  latitude  and 
longitude  covering  the  eastern  U.S.  (GALE  grid).  The  influence  of  two 
different  lateral  boundary  treatments,  namely  the  tendency  relaxation  scheme 
of  Perkey  and  Kreitzberg  (PK)  and  the  Davies  relaxation  scheme,  are  compared. 
For  integrations  with  initialized  fields,  the  NMC  analyzed  fields  interpolated 
to  the  model  grid  are  also  initialized  to  provide  balanced  boundary  values  and 
tendencies.  This  reduces  noise  in  the  boundary  zone,  which  can  propagate  into 
the  interior  of  the  domain  and  contaminate  our  test. 

In  the  static  initialization  procedure,  the  non-divergent  wind  is  first 
computed  for  the  analyzed  winds  on  the  pressure  surfaces,  by  solving  for  the 
streamfunction.  The  non-divergent  wind  and  the  analyzed  temperatures  are  then 
interpolated  to  the  sigma  coordinates  of  the  model.  A  diagnostic  relation  is 
then  derived  for  the  geopotential  on  the  sigma  surfaces  of  the  numerical 
model,  by  ignoring  the  tendency  of  divergence,  non-linear  vertical  advection 
and  friction.  The  initial  non-divergent  wind  and  analyzed  temperatures, 
interpolated  to  the  sigma  surfaces  are  used  to  compute  the  non-linear  forcing 
terms.  With  the  tangential  wind  defined  along  the  model  lateral  boundary  with 
the  C  grid,  consistent  boundary  conditions  for  the  normal  derivatives  of  the 
geopotential  are  easily  obtained  by  ignoring  the  tendencies  in  the  momentum 
equations . 

The  NRL  model  uses  the  split-explicit  scheme  to  integrate  in  time.  The 
scheme  is  compared  to  a  centered  difference  scheme  by  integrating  the  model  on 
the  GALE  grid.  The  split-explicit  scheme  has  been  shown  to  reduce  the 
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amplitude  of  the  unwanted  external  gravity  wave  oscillations  in  the  first 
three  to  four  hours  of  integration.  However,  a  typical  rapid  adjustment  with 
superimposed  oscillations  occurs  in  the  mid- troposphere  vertical  motion  in  the 
first  A  to  6  hours  of  integration.  The  new  static  initialization  procedure  is 
also  tested  on  the  GALE  grid.  By  interpolating  the  non-divergent  wind  and 
analyzed  temperature  to  the  model  sigma  surfaces,  the  amplitude  of  the  initial 
oscillations  of  the  surface  pressure  are  reduced  to  2  to  3  mb  and  are  largely 
damped  out  after  3  hours.  Using  the  non-divergent  wind  and  performing  a  static 
non-linear  balance  of  the  mass  field,  provides  a  balanced  initial  state, 
except  for  a  smooth  initial  adjustment  of  the  vertical  motion  in  the  first 
five  hours  or  less  of  integration  and  a  small  mean  drift  in  the  surface 
pressure.  For  the  varying  degrees  of  static  initialization,  similar  12  to  48 
hour  forecasts  are  produced  when  the  split-explicit  scheme  for  time 
integration  was  used. 

A  vertical  mode  initialization  scheme  following  that  of  Bourke  and 
McGregor  (1983)  has  been  developed  for  use  with  the  NRL  model.  Filtering 
conditions  are  applied  to  the  model  dynamical  equations  to  derive  to  linear 
diagnostic  equations  for  the  mass  divergence  and  geopotential,  which  are 
solved  iteratively  for  the  first  three  vertical  modes  of  the  numerical  model. 
These  modes  have  phase  speeds  which  are  much  faster  than  those  of 
meteorological  systems.  The  further  condition  that  the  linearized  potential 
vorticity  is  unchanged  by  the  procedure  is  required  to  compute  the  vorticity. 
The  observed  wind  and  temperature  is  first  interpolated  to  the  sigma  surfaces 
of  the  model.  The  iterative  procedure  is  then  used  to  compute  incremental 
changes  to  the  generalized  geopotential,  mass  divergence  and  vorticity  for  the 
first  three  vertical  modes  of  the  numerical  model.  As  is  customary,  we  keep 
the  geopotential,  temperature  and  pressure  fixed  at  the  lateral  boundaries  in 
the  scheme.  To  provide  a  boundary  condition  for  the  divergence  however,  an 
approximate  divergence  at  the  boundary  is  computed  using  the  thermodynamic 
equation.  In  our  scheme,  changes  in  the  tangential  wind  along  the  lateral 
boundaries  are  consistent  with  the  changes  in  the  vorticity  and  mass 
divergence  computed.  The  procedure  provides  a  balanced  vertical  motion  field 
and  produces  smaller  changes  to  the  initial  mass  and  wind  fields,  compared  to 
the  static  initialization.  The  scheme  is  tested  on  two  grids,  of  differing 
domain  size  and  grid  resolution.  Convergence  of  the  scheme  is  rapid  with  the 
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lower  resolution  US  grid,  with  three  iterations  of  the  scheme  being  sufficient 
for  convergence.  With  the  smaller  GALE  grid  of  higher  resolution  and  sloping 
topography  along  two  boundaries,  the  convergence  of  the  scheme  is  6lower  and 
in  fact  the  mass  divergence  didn’t  converge  for  the  third  mode.  However,  by 
smoothing  the  topography  in  a  boundary  zone  of  five  degrees,  the  convergence 
of  the  scheme  is  much  improved.  In  both  cases,  changes  in  the  mass  and  wind 
fields  are  still  small  after  three  iterations. 

Integrations  with  initial  conditions,  initialized  with  the  vertical  mode 
initialization  procedure,  prevent  gravity  wave  oscillations,  without  producing 
a  mean  drift  in  the  surface  pressure,  and  provide  a  balanced  vertical  motion 
field.  On  the  coarse  US  grid,  little  difference  is  found  between  integrations 
using  either  of  the  lateral  boundary  treatments.  However,  some  low  amplitude 
oscillations  in  the  surface  pressure  of  long  period  remain  in  the  interior  of 
the  domain  and  some  noise  is  generated  in  the  vertical  motion  at  the  lateral 
boundaries  when  the  Perkey  Kreitzberg  scheme  is  used.  On  the  smaller  GALE 
grid,  some  noise  is  produced  in  the  vertical  motion  in  the  boundary  zone  by 
both  schemes.  In  the  interior,  the  Davies  scheme  produces  a  smoother  variation 
of  the  surface  pressure.  When  the  Davies  scheme  is  used  with  integrations 
starting  from  uninitialized  data  and  boundary  values,  an  initial  shock  in  the 
surface  pressure  is  damped  in  the  first  four  hours.  However  a  small  drift  in 
the  surface  pressure  is  produced.  This  indicates  that  the  boundary  values  used 
with  the  Davies  scheme  should  be  as  balanced  as  possible  for  the  numerical 
model,  being  initialized  or  derived  from  integrations  on  a  larger  grid  or  with 
another  model.  Similar  12  to  48  hour  forecasts  are  again  produced  with  the 
various  experiments. 

For  grids  of  high  resolution  such  as  our  GALE  grid,  and  especially  when 
fine  scale  topography  is  used  along  the  boundary,  it  is  recommended  that  no 
more  than  three  iterations  of  the  vertical  mode  scheme  should  be  used  in 
practice.  For  grids  of  even  higher  resolution,  only  the  first  two  modes  may  be 
able  to  be  initialized,  with  possibly  no  more  than  three  iterations  used  with 
the  scheme.  Some  improvement  however,  can  be  expected  by  smoothing  the 
topography  in  the  boundary  zone  or  using  a  nested  model  to  provide  more 
accurate  boundary  values  for  the  mass  divergence  for  the  inner  nested  grids. 


48 


ACKNOWLEDGEMENTS 


The  first  author  was  supported  under  contract  number  N00014-86-C-2365  with 
NRL  and  the  second  author  was  supported  by  SPAWAR  and  NRL’s  basic  research 
program.  Preliminary  results  were  described  by  Sashegyi  et  al.  (1987)  in  the 
Proceedings  of  the  IAMAP  Symposium  on  Mesoscale  Analysis  and  Forecasting,  held 
17-19  August  1987  in  Vancouver,  Canada  and  by  Sashegyi  and  Madala  (1988)  in 
the  preprints  volume  of  the  Eight  Conference  on  Numerical  Weather  Prediction, 
which  was  held  22-26  February  1988  in  Baltimore,  Maryland. 


49 


APPENDIX :  SPLIT-EXPLICIT  SCHEME 


The  model  equations  in  matrix  notation  can  be  written 
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where  the  variables  are  defined  in  section  2.1,  and  the  non-linear,  Coriolis, 
friction  and  diabatic  terms  are  included  in  the  terms  on  the  right  hand  sides 
of  the  equations.  Integrating  (Al),  (A2),  (A3),  (AA)  and  (A5)  over  a  time  step 
2  At  we  obtain 
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where  the 
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and  p1  and  represent  grid  point  averages  as  defined  in  section  2.1.  The 

non-linear  and  Coriolis  terms  on  the  right  hand  sides  are  slowly  varying  so 
that  Ay  -  Au(t),  Ay  ■  Av(t),  A^  *  Aj(t).  For  an  explicit  time  step  choosing 
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#  2  #(t),  D  -  D(t),  gives  the  conventional  centered  difference  scheme  for 
time  integration.  We  use  this  explicit  difference  scheme  to  compute  first 


estimates  of  the  variables  at  time  t+At, 

giving 

pxuex(t+At)  -  p^uU-At)  +  2At  6  4(t) 

6  8  X 

3K 

2At  A  (t) 
u 

(A13) 

prvex(t+At)  -  p^v(t-At)  +  2At  6  #(t) 
s  s  y 

X 

2At  A  (t) 
v 

(A14) 

pgTeX(t+At)  -  pgT(t-At)  +  2At  M2  D(t) 

e 

2At  AT(t) 

(A15) 

PgX(t+At)  -  ps(t-At)  +  2At  N2T  D(t) 

ee 

0 

(A16) 

p  qex(t+At)  -  p  q(t-At) 

8  S 

- 

2At  Aq(t) 

(A17) 

Subtracting  these  equations  from  Eqs.  (A7),  (A8),  (A9),  (A10)  and  (All),  we 

can  then  write  for  the  corrected  variables 
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where  the  terms  on  the  right-hand  side  are  the  explicit  computations  of  the 
variables.  We  only  need  to  compute  the  corrections  ♦  -  f(t),  D  -  D(t)  for 
those  vertical  modes  whose  phase  speeds  relative  to  the  ground  U  +  i Xj[ 
are  greater  than  Ax/ (2At) ,  where  U  is  the  maximum  background  flow  speed.  It 
can  also  be  noted  that  the  specific  humidity  q  need  not  be  corrected.  To 
obtain  the  corrections  for  these  modes,  the  amplitudes  of  the  deviations 
♦  -  #(t),  D  -  D(t)  are  integrated  at  smaller  time  steps,  over  the  interval 
of  twice  the  large  time  step. 

To  integrate  the  deviations  we  require  equations  for  the  tendencies  of 
mass  divergence  and  generalized  geopotential.  By  taking  Sx  of  hy  times  Eq. 
(Al)  and  adding  6j  of  hx  times  Eq.  (A2)  we  obtain  the  equation  for  the 
divergence.  By  taking  times  Eq.  (A3)  and  adding  (RT*  +  f*)  times  Eq.  (A 4) 
an  equation  for  the  generalized  geopotential  is  obtained.  That  is 
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and  the  two  dimensional  Laplacian  is  defined  by  Eq.  (2.13)  in  section  2.1. 
In  terms  of  our  vertical  modes,  the  amplitudes  of  the  mass  divergence  and 
generalized  geopotential  are  d  *=  E_1D,  e  =  £“*♦.  The  amplitudes  of  the 
deviations  of  the  divergence  dkn  -  dk(t),  and  the  generalized  geopotential 
ekn  "  are  integrated  with  smaller  time  steps  given  by  Ark  =  At/m^. 

For  the  ten  layer  model,  the  deviations  are  computed  for  the  first  three  modes 
with  1%  given  by  8,  4  and  2  respectively.  Now  integrating  over  smaller  time 
steps  we  have 
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For  our  explicit  time  step,  applying  the  centered  difference  scheme  with  time 
step  2  At  to  Eqs.  (A23)  and  (A24),  or  from  Eqs.  (A13),  (A14),  (A15),  (A16) 
directly,  we  obtain 
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By  subtraction  we  obtain  the  equations  for  integrating  the  deviations 
dkn  -  dk(t),  ekn  -  ek(t)  , 
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For  the  integration  of  the  deviations  of  the  divergence,  a  lateral 
boundary  condition  is  required  for  the  generalized  geopotential.  The  boundary 
value  for  e^  -  ej^t)  is  computed  by  linear  interpolation  from  the  values  at 
t  -At  and  t.  Further  pragmatic  boundary  conditions  are  provided  by  reducing 
the  amplitude  and  phase  speed  of  the  deviations  of  the  divergence  and 
generalized  geopotential  in  the  boundary  zone  by  the  factor  1-a,  where  a  is 
defined  in  section  2.3,  for  the  two  different  lateral  boundary  formulations. 
That  is  to  reduce  the  phase  speed  c^  in  the  boundary  zone,  we  multiply  in 
Eq.  (A31)  and  the  term  in  Eq.(A30)  by  (1-a).  The  amplitudes  are  reduced  in 
the  boundary  zone  by  multiplying  the  correction  terms  in  Eqs.  (A20)  and  (A21) 
by  (1-a)  and  those  in  Eqs  (A18)  and  (A19)  by  (l-au)  and  (l-av),  respectively. 


where  (t+At)  = 
and  ekX(t+At)  = 
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SIGMA 


FIG.  1.  The  vertical  structure  of  the  modes  versus  sigma  level  in  a  ten 
layer  model  for  the  case  of  a  basic  state  at  rest  with  the  mean  temperature 
profile  defined  in  Table  1. 
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FIG.  2.  The  horizontal  grid  stencil  for  the  Arakawa  C  grid  as  applied  to 
the  NRL  model,  showing  the  lateral  boundary  grid  points,  u  and  v  represent  the 
horizontal  wind  components.  The  surface  pressure  ps,  the  temperature  T,  the 
geopotential  specific  humidity  q,  the  mass  divergence  D  and  vertical 

motion  6  are  defined  at  the  same  horizontal  grid  point  as  the  generalized 
geopotential  ♦. 
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FIG.  U.  (a)  The  analyzed  sea-level  pressure,  wir.d  and  temperature  at  1000 
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FIG.  5.  The  horizontal  grid  stencil  used  to  compute  the  changes  in  the 
wind  field,  given  the  changes  in  the  mass  divergence  and  vorticity.  The  stream 
function  change  At  is  defined  in  the  interior  at  the  same  grid  point  as  the 
vorticity  while  the  velocity  potential  change  Aj  is  defined  at  the  same 
grid  point  as  the  mass  divergence  D.  The  lateral  boundary  conditions  for 
and  At  are  indicated. 
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FIG.  6.  The  smoothed  model  topography  for  (a)  the  US  grid,  (b)  for  the 
GALE  grid  and  (c)  for  the  GALE  grid  with  the  topography  further  smoothed  in  a 
five  degree  zone  around  the  lateral  boundary.  The  horizontal  resolution  is  2° 
longitude  by  1.5°  latitude  for  the  US  grid  and  0.5°  in  latitude  and  longitude 
for  the  GALE  grid.  The  contour  intervals  are  500  m  for  elevations  above  1000 
m,  and  250  m  for  those  below  1000  m. 
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FIG.  7.  The  root-mean- square  (rms)  changes  in  the  amplitudes  on  the  US 
grid  of  (a)  the  mass  divergence  Ad  in  dynes  cm-2  s”1,  (b)  the  generalized 
geopotential  Ae  in  units  of  10^  dynes  s-2,  and  (c)  the  vorticity  Av  in  dynes 
cm”2  s”1,  versus  the  iteration  number  in  the  vertical  mode  initialization 
scheme,  for  each  of  the  first  three  vertical  modes.  In  (d)  is  shown  the  rms 
change  Aps  in  the  surface  pressure  in  mb  versus  the  iteration  number.  The  rms 
values  are  averaged  for  all  the  analyses  of  the  week  of  12Z  January  23  to  12Z 
January  29,  1986.  (e)  The  rms  change  in  the  amplitude  of  the  mass  divergence 
Ad  in  dynes  cm-2  s”^-  versus  iteration  number  for  the  highest  mode  initialized, 
as  the  number  of  modes  initialized  is  varied  from  one  to  six  modes  for  the 
analysis  on  12Z  January  23,  1986  only. 
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FIG.  7.  The  root-mean-square  (rms)  changes  in  the  amplitudes  on  the  US 
grid  of  (a)  the  mass  divergence  Ad  in  dynes  cm-2  s'^,  (b)  the  generalized 
geopotential  Ae  in  units  of  1011  dynes  s"2,  and  (c)  the  vorticity  Lv  in  dynes 
cm*2  s’1,  versus  the  iteration  number  in  the  vertical  mode  initialization 
scheme,  for  each  of  the  first  three  vertical  modes.  In  (d)  is  shown  the  rms 
change  Aps  in  the  surface  pressure  in  mb  versus  the  iteration  number.  The  rms 
values  are  averaged  for  all  the  analyses  of  the  week  of  12Z  January  23  to  12Z 
January  29,  1986.  (e)  The  rms  change  in  the  amplitude  of  the  mass  divergence 
Ad  in  dynes  cm-2  s_1  versus  iteration  number  for  the  highest  mode  initialized, 
as  the  number  of  modes  initialized  is  varied  from  one  to  six  modes  for  the 
analysis  on  12Z  January  23,  1986  only. 
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.ITERATION  NUMBER 


(C) 


FIG.  12.  The  root-mean-square  (rms)  changes  on  the  GALE  grid  in  the 
amplitudes  of  (a)  the  mass  divergence,  and  (b)  the  vorticity  versus  iteration 
number  in  the  vertical  mode  initialization  scheme  for  each  of  the  first  three 
vertical  modes,  averaged  for  all  the  analyses  for  the  week  of  12Z  January  23 
to  12Z  January  29,  1986.  (c)  Comparing  the  rms  changes  in  the  amplitude  of  the 
mass  divergence  when  the  GALE  topography  (see  Fig.  6b)  is  used  (solid  lines) 
to  that  when  the  GALE  topography,  which  has  been  further  smoothed  (see  Fig. 
6c)  in  a  five  degree  zone  around  the  lateral  boundary,  is  used  (dashed  lines). 
The  analyses  only  on  12Z  January  23  (labelled  1)  and  12Z  January  24,  1986 
(labelled  II)  are  used. 
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FIG.  13.  The  analyzed  (a)  sea-level  pressure  and  (b)  wind  field  on  the 
GALE  grid  at  the  model  sigma  level  a  =  0.25  for  12Z  January  23,  1986.  The 
contours  are  as  in  Fig.  8a  and  b.  The  maximum  wind  vector  shown  is  54  m  s'1, 
(c)  The  surface  pressure  change  in  mb,  after  the  vertical  mode  initialization. 
Contours  of  surface  pressure  change  are  every  0.25  mb. 


SURFACE  PRESSURE  AT  90W.  35N  (b)  VERTICAL  M0T10N  AT  LEVEL  5  AND  901-1 .  35N 
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FIG.  15.  Time  series  of  (a)  surface  pressure  in  mb  and  (b)  the  vertical 
motion  (in  sigma  coordinates)  in  hour"1  at  sigma  level  O  =  0.5,  at  a  grid 
point  on  the  GALE  grid  at  90°W  and  35°N,  during  the  first  12  hours  of 
integration.  Curves  are  for  uninitialized  initial  conditions  for  an  explicit 
integration  in  time  (Curve  A)  and  a  split-explicit  integration  (Curve  B). 


SURFACE  PRESSURE  AT  90W,  35N  VERT1CRL  M0T10N  AT  LEVEL  5  AND  90U.  35N 
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point  on  the  GALE  grid  at  90°W  and  35°N,  during  the  first  12  hours  of 
integration.  Curves  are  for  uninitialized  initial  conditions  for  an  explicit 
integration  in  time  (Curve  A);  for  a  split-explicit  integration  with  a  non- 
divergent  initial  condition  (Curve  C)  and  for  split-explicit  integrations  with 
a  non-linear  mass  balance  (Curve  D) . 
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(a)  SURFACE  PRESSURE  AT  90W.  35N  (b)  SURFACE  PRESSURE  AT  90H,  35N 
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VERTICAL  H0T10N  AT  LEVEL  5  AND  90W.  35N  (d)  VERTICAL  M0TI0N  AT  LEVEL  5  AND  90W,  35N 


VERTICAL  M0TI0N  AT  LEVEL  5  AND  46H.  40N  (d)  VERTICAL  M0TI0N  AT  LEVEL  5  ANO  4QN 
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VERTICAL  M0TI0N  AT  LEVEL  5  ANO  90W.  35N 
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SURFACE  PRESSURE  AT  83W.  35N 
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VERTICAL  M0TI0N  AT  LEVEL  5  AND  90W,  35N  (d)  VERTICAL  M0TI0N  AT  LEVEL  5  AND  83W.  35N 


hours  of  integration,  (a)  shows  the  surface  pressure  in  mb  at  90  W  and  3b  N 
and  (b)  at  and  35°N.  (c)  shows  the  vertical  motion  in  hour"1  at  90°W  and 
35°N  and  (b)  at  83°W  and  35°N.  Curves  are  for  uninitialized  initial  conditions 
(Curve  D)  and  for  initialized  initial  conditions  (Curve  E),  using  the  Davies 
lateral  boundary  scheme. 


SURFACE  PRESSURE  AT  59W.  40N  W  SURFACE  PRESSURE  AT  58W.  40N 
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at  59°W  and  40°N  and  (b)  at  58°W  and  4CTN.  (c)  snows  tne  vertical  motion  in 
hour'1  at  59°W  and  40°N  and  (d)  at  58°W  and  40°N.  Curves  are  for  uninitialized 
initial  conditions  (Curve  B)  and  for  initialized  initial  conditions  (Curve  C), 
using  the  Perkey  Kreitzberg  lateral  boundary  scheme.  Curve  E  is  for 
initialized  initial  conditions  with  the  Davies  lateral  boundary  scheme. 


59N,  'ION  VERTICAL  M0TI0N  AT  LEVEL  5  ANO  50M.  'ION 


SURFACE  PRESSURE  AT  59W.  40N  (b)  SURFACE  PRESSURE  AT  58W,  AON 


at  59°W  and  AO°N  and  (b)  at  58°W  and  AO°N.  (c)  shows  the  vertical  motion  in 
hour'1  at  59°W  and  AO°N  and  (d)  at  58°W  and  40°N.  Curves  are  for  uninitialized 
initial  conditions  (Curve  D)  and  for  initialized  initial  conditions  (Curve  E), 
using  the  Davies  lateral  boundary  scheme. 


VERTICAL  M0TI0N  AT  LEVEL  5  ANO  59W.  40N  (d)  VERTICAL  H0TI0N  AT  LEVEL  5  AND  58W.  ION 
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